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ABSTRACT 

The  low  frequency  scattering  of  electromagnetic  and  acoustic  waves  by 
rotationaliy  symmetric  bodies  is  considered.  By  concentrating  on  certain 
quantities  such  as  the  normalised  component  of  the  induced  electric  and  mag¬ 
netic  dipole  moments,  it  is  shown  how  the  first  one  or  two  terms  in  the  far 
zone  scattered  fields  can  be  expressed  In  terms  of  quantities  which  are  func¬ 
tions  only  of  the  geometry  of  the  body.  Each  of  these  is  the  weighted  integral 
of  an  elementary  potential  function  which  can  be  found  by  solving  an  integral 
equation.  A  computer  program  has  been  written  to  solve  the  appropriate 
equations  by  the  moment  method,  and  for  calculating  the  dipole  moments,  the 
electrostatic  capacity,  and  a  further  quantity  related  to  the  capacity.  The 
program  is  described  and  related  data  are  presented. 
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i.  INTRODUCTION 

When  a  plane  electromagnetic  wave  is  incident  on  a  finite  perfectly  con¬ 
ducting  body,  or  a  plane  acoustic  wave  incident  on  a  finite  acoustically  soft  or 
hard  body,  die  scattered  field  in  the  far  zone  can  be  expanded  in  a  power 
series  in  the  wave  number  k  if  k  is  sufficiently  sm«ll.  The  determination  of 
the  first  few  terms  in  these  series  requires  the  solution  of  certain  elementary 
potential  problems.  We  here  consider  the  potential  problems  associated  with 
the  first  ( Rayleigh )  term  in  the  electromagnetic  expansion  and  the  first  two 
terms  in  each  of  the  acoustic  expansions,  and  show  how  in  the  case  of  a  singly 
connected  body  of  revolution  all  of  these  terms  can  be  deduced  from  the  solu¬ 
tions  of  just  five  potential  problems.  If  the  body  is  not  singly  connected,  only 
the  axial  component  of  the  induced  electric  dipole  moment  is  affected,  and  for 
a  body  consisting  of  two  separate  parts,  an  expression  for  the  modified  compo¬ 
nent  is  obtained. 

Each  potential  satisfies  a  simple  integral  equation.  Computer  programs 
are  described  for  solving  the  equations  by  the  moment  method,  and  since  most 
of  the  equations  are  of  first  order  type,  the  computational  procedures  are  rather 
similar  to  those  of  Mautz  and  Harrington  ( 1970).  The  appropriate  elements  of 
the  electric  and  magnetic  polarisability  tensors  (Keller  et  al,  1972)  are  then 
computed,  along  with  the  electrostatic  capacity  and  a  quantity  7  related  to  this, 
and  these  are  sufficient  to  specify  the  electromagnetic  and  acoustic  scattering 
for  any  direction  of  plane  wave  incidence  and  any  direction  of  scattering.  For 
relatively  simple  geometries,  tee  entire  computation  takes  about  3  seconds  on  an 
IBM  360  computer. 

In  our  presentation  we  first  examine  (Section  2)  the  problem  of  a  plane 
electromagnetic  wave  of  arbitrary  polarisation  and  incidence  direction,  and 
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isolate  the  potentials  neoessary  for  a  complete  description  of  the  leading  term 
in  the  far  zone  scattered  field.  This  is  followed  (Section  3)  by  similar  treat" 
ments  of  the  acoustic  problems,  but  here  we  seek  the  firsc  two  terms  in  the 
expansions.  In  Section  4  the  integral  equations  satisfied  by  the  potentials  are 
cast  into  forms  appropriate  to  digital  solution,  and  the  manner  in  which  the 
body  is  specified  is  also  described.  Section  5  is  concerned  with  various  aspects 
of  the  computer  program,  a  complete  listing  of  which  is  given  in  the  Appendix, 
and  some  of  the  numerical  results  obtained  so  far  are  presented  in  Section  6. 

This  work  was  supported  in  part  by  the  National  Science  Foundation  under 
Grant  GP-25321. 
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2.  PERFEC1 


2.  i  FORMULATION 

I/et  B  be  a  finite,  closed,  perfectly  conducting  body  of  revolution  about 
‘he  z  axU  of  a  rectangular  Cartesian  coordinate  system  (x,y,  z).  In  terms  of 
C .  j  cylindrical  polar  coordinates  (p,0,  z)  where 


f  2  2 

p  v/'*  +y 


0  *  arctan  ^  , 


the  surface  will  be  described  by  the  equation 
P  *p(z) 

where  p  can  be  a  multivalued  function  of  z  as,  for  example,  in  the  case  of 
a  disk  or  a  re-entrant  shape,  but  is  never  infinite  and  is  zero  outside  some 
interval  in  z .  Let  £  be  the  radius  vector  to  an  arbitrary  T>oint  in  the  domain 
V  exterior  to  B  and  let  n  be  a  unit  vector  normal  to  the  surface  drawn  into  V. 

A  linearly  polarised  electromagnetic  wave  is  incident  with  electric  and 
magnetic  vectors 

A 

J  ik  K  r 

E*  a  e 


„i  ..  A  i  k  k,  r 
H  *  Ybe 


A  A  A  A  A  A 

where  k,  a  and  b  are  mutually  perpendicular  unit  vectors  such  that  b*-  kAn  ; 

Y  is  the  intrinsic  admittance  of  the  homogeneous  isotropic  medium  ( of  permit- 

**iu( 

tivity  e)  exterior  to  B  and  a  time  factor  e  has  been  suppressed. 
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B  & 

For  k  small  but  kr  largo,  the  resulting  scattered  field  E  ,  H  can  be 
written  as  ( Kleinman,  1965 ) 


_s  .  .2  r 

—  ^  4 » r 


A  ,  A  . 

E  rA<rA£>  + 


K 


ikr 

e 

4x  r 


2  f  Y 

1  l7 


(?AE)-fA<fA 


y 

m)J 


(2) 


where  £  and  jn  are  the  electric  and  magnetic  dipole  moments  respectively. 
As  shown  by  Keller  et  al  ( 1972 ) , 


(P 


**«  <Pa?+  (P33 


m  -  -y  <Mub  f  (M33 


“  Pjj)  (tL.  Z)  Z  ^ 

.  ,  A  A  A  \ 
-  )  (b.  z)  z  J 


(3) 


whe  i  P^,  and  are  functions  only  of  the  geometry  of  the  body. 

For  a  given  body,  P^,  P33,  and  M33  are  conrtants  whose  values  are 
as  follows: 


(i) 

Pn’  JJ  x 

(4) 

where 

B 

■$ ,  is  such  that 

<1 

fO 

y* 

H 

O 

in  V 

!i  ■ * 

on  B 

(5) 
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(iv) 


|(r“2) 


//s  5 


as  r  —400  . 


n.  z  (z  dS 


B 


(in 


whore  *5**  is  such  that 
o 


V2  4,  *  0  In  V 


!±3_  „  li 
9n  Bn 


on  B 


0(r"2) 


as  r  -400 . 


02) 


Although  the  values  assumed  by  the  potential  function  on  B  are  quite  distinct 

T  a$  ad  6 

from  those  of  •y.  ,  J X  and  *3  t  nevertheless,  as  shown  by  Karp 

1  an  an 

(1956)  and  Payne  ( 1956), 


2  P11 


( 13) 


This  obviates  the  need  for  solving  the  potential  problem  ( iv)  if  the  only  purpose 
for  finding  ijr3  is  to  calculate  M33  . 

There  is  one  other  electromagnetic  quantity  of  interest  and  this  is  the 
electrostatic  capacity  C  of  the  body  In  isolation.  If  die  body  is  raised  to  die 
potential  unity,  the  surface  charge  density  is 
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where  $c  is  an  exterior  potential  function  satisfying  the  boundary  condition 
§L  *  1  onB  .  ( 


The  electrostatic  capacity  is  then  equal  to  the  total  charge  induced  on  the  sur¬ 


face  and  is 


C  *  -e 


a 


Note,  however,  that  if  all  portions  of  the  surface  are  not  in  electrical  contact, 
with  one  another,  charge  can  no  longer  flow  freely  over  the  entire  surface,  and 
additional  ( mutual )  capacities  can  be  defined,  in  particular,  such  electrical 
separation  has  a  profound  effect  cm  the  calculation  of  P  ,  and  the  modifications 
that  result  when  the  surface  is  disjoint  are  discussed  in  Section  2. 

The  five  quantities  listed  in  sqs.  ( 4 ),  ( S ),  ( 9  ),  ( 11 ),  and  ( 16 )  can  be 
computed  by  solving  five  separate  potential  problems  of  a  rather  standard  nature, 
and  the  manner  in  which  this  is  done  is  as  follows. 

2 

Let  V  be  some  potential  function  satisfying  V  V  »  0  outside  and  on  B, 

3  8 

and  let  V  be  the  regular  part  of  V.  V  is  therefore  an  exterior  potential 


and  we  can  regard 


Y  -  V8  *  V* 


as  an  incident  potential.  Green's  theorem  applied  to  the  function  V  in  the  region 
V  then  yields 


V  (g)  *  V*< 


-  jr f?vt*',}dS 


7 


<4YJ  *1*0 , 
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where  R*|r  -  r’J  . 

If  the  boundary  condition  on  the  potential  V  is 
V(r)  «  0  ,  r  on  B, 

eq.  ( 18 )  reduces  to 

VCD  -  v‘(X)  -  fj±  -g;,  VtrldS1  . 

B 


( 19) 


(20) 


a ‘id  since  the  integral  exists  for  all  r  Including;  points  on  B,  we  can  allow  jr 
to  lie  on  B  and  apply  the  boundary  condition  ( 19 )  to  obtain 


V*(  r) 


-h  jjkh  V(I* 


)  ds'  , 


B 


whica  is  an  integral  equation  of  the  first  kind  for 


9  V 
9n 


If,  on  the  other  hand,  the  boundary  condition  on  the  potential  V  is 


(21) 


bn 


V(r)  ■  0  ,  r  on  B, 


( 22 ) 


eq.  ( 18 )  reduces  to 


V( 


£)*v‘(r)+i-  //Wj&GJm*  . 


B 


(23) 


and  because  of  the  non-integrable  singularity  of  ~r  at  r*  *  £  , 

eq,  ( 23 )  is  valid  as  it  stands  only  if  r  is  not  on  B.  To  obtain  an  integral 
equation  for  V,  we  allow  £  to  approach  a  point  on  B  in  the  direction  of  the 


1  . 


i 


t  «e 

“  £  * 
s  s 


6  fe 


I  * 
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inward  normal,  in  which  case  it  can  be  sho.vn  that 


lim 

r-*B 


//  v(i,)  h  (i)dS'-  2'rV<£>  +//v<£  >  i^(i)ds' 


where  the  bar  across  the  integral  signs  denotes  the  Cauchy  principal  value.  Hence 

3  /l  \  c» 

I  -  3?S-US 

B 


V(r)  »  2Vi(r}+  ~ 


lR/ 


dS*  (24) 


for  r  on  B,  which  is  an  integral  equation  of  the  second  kind  for  V. 


2. 2  PROCEDURE  FOR  P 


11 


The  solution  of  problem  (i)  and,  hence,  the  computation  of  F  is  a 
straightforward  application  of  the  integral  equation  ( 21 ) . 

If 


ti  X 

«  x  ,  V  * 


then 


V  *  x 


-  $  1  *  v 1  ( say ) 


(25) 


with  V  »  0  on  B.  Since  x  ■  p  cos  0  and  the  equation  of  the  surface  is  indepen¬ 
dent  of  0,  the  potential  must  everywhere  have  the  same  0  dependence 
as  V  *,  implying 


V^r)  «  Vj(p,  z)  cos  0  . 
dVl 

This  ie  true  also  of  ~ —  ,  and  we  can  therefore  write 

on 


(26) 


^ggg 
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3n 


»  v1(£*)  *  v  1(p'  .  *' )  coa  0* 


«  T1{s,>coe0' 


where  s'  is  arc  length  along  a  profile  cf  the  body.  Moreover, 


(2?) 


with 


R  *  ju-p' 


^  *  0  “  0‘ 


2  9  V/2 

)  +  (  z  -  z  * )  +  2  p  p  *  ( 1  -  coe  )  i  <  28 ) 


and  since 


dS'  *  p’ d  0*ds*  . 

the  integral  equation  ( 21  /  now  takes  the  form 

2a  js 


p  cos  0 


■7 rf  /  ^  P’^’ds' 

•D  *0 


2  a  ts 


17  |  /  T1  {  8  *  >  ^^p’Wds' 

»0  *0 


8 

1__  f  p’XjTj  (s’)ds' cos  0 


where  the  kerrel  is 


K 


a 

l  *  Ki<p‘  zi  *>*•  -*>  *  J 

•a 


coat'/ 


d#  ( 29 ) 
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and  the  integrrtion  is  along  the  profile  of  the  body.  The  integral  equation 
for  T^Cs)  is  therefore 


s 

/. 


p  f  T^(s  * )  ds  '  *  2  ?p 


(30) 


which  can  be  solved  to  date  r  mine  T  ( s )  .  In  terms  of  this  quantity 


11 


// 

B 


x  -3TV(r)dS 


2r  s 


*  /  / 
* 0  «*0 


p2  cos2  0T(s)d0ds 


which  reduces  to 


11 


t  2 

T  I  o  T,(s)ds 

•an 


(31) 


2. 3  PROCEDURE  FOR  P  AND  C 

The  solution  of  problem  ( ii ).  leading  to  the  calculation  cf  P^,  involves 
two  successive  applications  of  the  integral  equation  (21 ) . 

In  the  first  case  we  consider  an  incident  potential 

V1  -  V*  «  X  (32) 

and  seek  the  corresponding  total  potential  V  satisfying  the  boundary  condition 


11 


013630-9-T 


,i 


V  *  0  on  B.  Since  V  is  everywhere  independent  of  0  ,  it  follows  that 
2  3V2 

V2  and  g  ■  ara  likewise  0  independent.  We  can  therefore  write 


av, 


IT-  *  V*'> 


and  die  integral  equation  ( 21 )  now  becomes 


2  jr  s 


« /  / 

•0  *0 


T2(8,)  R  P’d^d8’ 


T7  f  P'  K„T,(.')ds> 

•n 


2x  I  “0 

*0 


where  die  kernel  is 


Kq  2  Kq  ( p,  z ;  p 1 ,  z  *) 


/ 


d 

R 


The  integral  equation  from  which  to  determine  T^(  s )  is  therefore 


8 

l 


p’  KgT^s'Jds'  -  2sz  . 
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(33) 


(?4: 


(35) 
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The  second  of  the  two  basic  problems  is  that  in  which  the  incident  potential 


is 


V1  «  vl  -  1  . 


(36) 


We  again  seek  the  total  potential  satisfying  the  boundary  condition  V  *  0 
on  B,  and  writing 


3V 

IT  *  V*  >  • 


(37) 


the  integral  equation  ( 2.1 )  takes  the  form 


s 

I, 


p,K0Tg(s*)  ds*  «  2*  , 


(38) 


from  which  T  ( a )  can  be  found. 

J 

In  problem  ( ii ) ,  however. 


V1* 


z  +  7  and  V8  «  - 


Thus. 


V1  .  V*  +  7V3‘ 


(39) 


implying 


V  "  V2  +  7V3  * 


(40) 
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and  if  we  write 


«  T(b')  , 


T(s)  *  T2<  8  >  +  YT3(«} 


where  T^s)  and  T^Cs)  are  the  solutions  of  the  Integral  equations  (35)  and 
( 38 )  respectively.  The  constant  y  is  determined  by  the  condition  ( 8 )  for  zero 
total  induced  charge  on  B,  viz.  * 


// 


( v1  - V)  dS  «  0  . 


//inrds  ‘  Hi i  <i+*>dS-°  . 


as  can  be  seen  by  application  of  the  divergence  theorem;  moreover 


Jp  “  i  l 


T(  s)  pdjf  dS 


f 

«  2a  J  pT(s) 


*  We  are  here  assuming  that  the  surface  is  not  disjoint. 
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s  8 

»  2rJ^  p  T2(s)d8  +  7^  p  Tgfsjds 


and  hence,  by  virtue  of  eqs.  ( 43 )  and  ( 44 ), 


p  T2<  8 )  d  8 


p  Tg  (  8 )  ds 


Since  Tg(  b)  and  Tg( a )  can  be  found  from  the  integral  equations  ( 35 ) 
and  ( 38 ) ,  the  constant  7  given  in  eq.  ( 45 )  now  completes  the  specification 
of  the  surface  field  T(  a  ),  and  in  terms  of  T(  s ) 


■// 


1  -5T  v(£,ds 


*  2a  I  spT(s)  ds  . 


Hence 


8 

-  2a  | 

Jo 


zp  T2<s)  ds  +  2  a  7  I  zpTg{s)ds  . 


A  valuable  by-product  of  the  above  analysis  is  the  electrostatic  capacity 
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C  defined  in  eq.  (16 ).  This  fact  is  apparent  on  recalling  that  the  determination 
of  C  requires  us  to  find  the  exterior  potential  satisfying  the  boundary 
condition  (  15 )  on  B,  and  this  can  be  accomplished  using  the  integral  equation 
(21)  with 
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C  /(2  *  e),  will  oh  ensures  that  1  can  never  be  infinite. 

6 

Some  simplification  of  the  preceding  results  is  possible.  Since  V  and 

8 

V  are  both  exterior  potentials.  Green’s  theorem  can  be  applied  to  the  region 

V  exterior  to  B  to  yield  the  reciprocity  relation 

6  8V3  t  [  s  SVt 

V 2  ~iras  ’  JJV !  1T<S 

B  B 


( Van  Bladel,  1958).  But 


and 


dn 


T2<*> 


A  A 

n  .  z 


9n 


T3(s)  . 


and  from  the  boundary  ccrditions  on 


and  V3 


o 


M  -1 


on  B.  Substituting  these  into  eq.  ( 51 ) ,  we  have 


zT3(s)dS  * 


B 

which  reduces  to 


8 


pT2(s)ds  « 


j l  f"T2(S^  ”n  * 


s 


zp  T3  (  s )  d  s  . 


(52) 
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With  the  aid  of  this  result,  the  constant  7  of  eq.  ( 45 )  can  be  expressed 
in  terms  of  the  surface  field  ( s )  alone: 


s 

/ 


zpT^(s)d  s 


7  *  “ 


(53) 


l 


p  T^( s ) ds 


i.e. 


s 

l 


T-  -  ^  |  zpT3(s)ds  , 


( 54) 


but  whilst  this  reduces  from  four  to  three  the  number  of  separate  surface  field 

integrations  involved  in  the  calculation  of  P^,  there  is  no  way  to  avoid _ 

entirely  the  determination  of  the  surface  field  T  ( s ) .  Indeed,  the  simplest 
expression  for  P  is 

uw 


s 

/ 


2  C 


P33  *  2*  I  T9(s)ds  -  7  7 


(55) 


2.4  PROCEDURE  FOR  M 


11 


The  solution  of  problem  ( iii)  leads  to  die  calculation  of  and  is  a 
straightforward  application  of  the  integral  equation  ( 24 ) .  If 


V  *  x 


and 


"  a  -%  ’ 
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with  ^  *  0  -  0‘  as  before.  Hence 


2  17 

l 


008  ^  &  (s' 


^  2ir 

jd  0  *  »  cos  0  ^ 


COS  i> 


'\ 


L 


an*  R 


d  (/; 


*  2  cos  0  jp  cos  a'  «2  +  j^(  z  ’  -  z )  sin  a 1  -  p  *  cos  a  ^  l  ^ 


where 


ni  *  °1 (p*  z;pf*  z'* 


2  2 


i 

i 


cos 

id 

3 

R 

2  , 

COS 

(61) 


(  G2  ) 


R 


The  integral  equation  from  which  to  determine  V4  ( s )  is  now 
V4(b)«  2p+“  £  v4(s')  fpcoea'fl2+  (te*  -  z)  sir. e* -p1  cos  oQ  n^jp*  ds! 


(63) 


1 
.  u 


i  3 
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i.e. 


P  V  {  s  )  COB  at  dfi 
4 


2.5  PROCEDURE  FOR  Mg3 


(64) 


Although  it  is  not  necessary  to  compute  M  oirectly  because  of  die 
relation  (13),  the  integral  equation  which  the  corresponding  potential  satis¬ 
fies  must  be  solved  if  the  second  term  in  the  low  frequency  expansion  lor  an 
acoustically  hard  body  is  to  be  evaluated.  It  is  therefore  appropriate  to  des¬ 
cribe  the  determination  of  this  potential  function  here. 

Once  again  we  have  a  straightforward  application  of  h.e  mtegnl  equation 
(24).  If 

V  1  *  2  and  V  3  * 

•  ? 

O 

then 

V  =»  z  -  Y  *  V„  ( say  /  ( 65) 

J  0 


9V 

with  5/dn  ■  0  on  B.  Since  V_  must  be  independent  cf  jt,  eq.  (24) 

5 

implies 


2  s-  8 


V5(S 


}  9  2Z+  **  fo£V 5<s,)  ' 


But 


2s 


2s 


l  *■'  (") d0,  *  I  *  w 


d 
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.  2^cob«'  ju'-sj.tas'-p'cosa  ]"•) 
where  f2  is  at  defined  in  eq.  ( 61 )  and 


ft0  *  V  P.sjp'.z'l 


l 


M. 

R3 


(  <>(>  ) 


The  integral  equation  from  which  to  determine  V,.  { s }  is  therefore 


V(a)*2z  +  — 
5  x 


and  we  note  in  passing  that 


0 

Vgte’J^pcoe  a'  flj  +  ^z*- 


z)  sino*  -  p  *  cos  a 


M33- 


// 


n.z  V  (e) dS 

5 


*]  nAp’ds* 

J  J  (67) 


B 


i.e. 


2*  s 

ll 

8 


sin  ( s )  p  d  0  d  s  , 


/ 


M33  *  “2a|pVc(s)3ino'ds 


(68) 


2.6  DISJOINT  SURFACES 

So  far  it  has  bean  assumed  that  all  portione  of  the  surface  Are  in  electrical 
contact  with  one  another,  and  if  this  requirement  is  not  met,  the  analysis  is  no 
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longer  valid.  Thus,  for  example,  at  application  of  the  above  procedures  to  a  body 
consisting  of  two  separate  spheres  leads  instead  to  the  solution  for  the  two  spheres 
joined  by  an  infinitesimal  wire  along  the  axis  of  symmetry,  and  though  the  presence 
of  the  wire  (producing  elsetrical  contact)  does  not  affect  the  values  of  and 
(  and  hence  M,^ ,  by  virtue  of  eq.  13 ),  it  does  have  a  profound  effect  on 
P  .  This  is  not  unnatural  since  P  is  proportional  to  the  longitudinal  ( z ) 

OO  00 

component  of  the  induceo  electric  dipole  moment. 

The  breakdown  in  our  formulation  when  B  has  several  distinct  parts  stems 
from  the  imposition  of  the  zero  induced  charge  c rite; ion  (9).  If  charge  cannot 
flow  freely  between  the  n  parts  B  ,  B  f _ ,  B  ,  eq.  (8)  must  be  replaced 

A  I  n 

by  the  n  equations 


i  *  1,  2, - -  n.  { 69 ) 


Since  this  obviously  alfects  only  the  potential  §  ^  and  leaves  the  procedure 
(and  results)  for  P, ,,  M,,  and  M__  unchanged,  our  efforts  will  be  directed 

ll  11  OO 

at  P  alone  with  the  objective  of  finding  an  approach  which  is  apr.'josble  when 

OO 

B  consists  of  just  two  electrically  Isolated  portions  B^  and  B2  .  So  that  we  may 
use  to  the  fullest  extent  the  work  that  we  have  already  done,  it  is  desirable  to 
have  this  new  approach  as  similar  as  possible  to  that  appropriate  when  the  two 
portions  are  electrically  connected. 

By  analogy  with  problem  (ii)  of  Section  2. 1,  the  task  is  to  find  an  exterior 
potential  satisfying  the  equation  V  “  ^  »  0  in  the  domain  V  exterior 

O  O 

to  B ,  together  with  the  boundary  conditions 


♦  3  '  1  +  Y1 

on  B, 

1 

(70) 

$3  •  i  +  f2 

on  B 

2 

(71) 
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where  the  constants  7,  aod  7^  are  such  that 


(72) 


(73) 


The  quantity  P  ifi  then  given  by  eq.  ( 6  )  as  before. 

OO 

Because  the  boundary  conditions  on  B  and  B„  differ,  it  ie  no  longer 
convenient  to  think  in  ter  ns  of  incident  and  total  potentials,  with  the  difference 
representing  the  desired  exterior  potential.  Let  uc  therefore  consider  the  basic- 
potential  problem  in  which  \  ^  ^  is  an  exterior  potential  satisfying  the  boundary 
condition 


(1) 


on  Bj 
on  B2 


(74) 


By  application  of  Green's  theorem  to  the  domain  \) ,  we  have 

*,"v> -t  *  h 4- 

B  ^ 


tffMih  -hjjihi 

B.  B 


3a)(r')dS' 
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and  the  first  integral  Is  identically  zero  since  it  itself  a  closed  surface. 
If,  now,  r  is  allowed  to  approach  B ,  application  of  the  boundary  condition 
(74)  gives 


77  If  n  h 

B  ,  r  on  B 


("5) 


which  is  an  integral  equation  from  which  to  determine 


3 


(1) 


>3b<  ■  itca° 


x  (  1 )  3  6 

be  simplified  somewhat  by  observing  that  „  and.  hence,  T  3 


(1) 


are 


3n 


are  independent  of  the  azimuthal  coordinate  0 .  When  the  0  integration  is 
performed,  eq.  ( 75)  reduces  to 


p‘  K0T3(1)(s,)ds’ 


j  2  ,  r  on 

^0  ,  r  on  B9 


c.  f.  eq.  (3$),  where 


(76) 


(77) 


and  Kq  Is  the  kernel  defined  in  eq.  ( 34 ) .  It  will  be  noted  that  the  integration 

in  ( 7fi )  is  over  the  entire  profile  of  the  body  B  *  B  *  B„  . 

(  2 )  12 
Similarly,  if  $  3  is  an  exterior  potential  function  satisfying  the 

boundary  condition 


$ 


(2) 

3 


fo  on  Bj 
[l  on  B2 


(78) 
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then 


f  P'K0T3(2>U')dS' 

Jn 


*  a— 


on  B . 


2t  .  r  on  B„ 


where 


t3<2><»,> 


X,  *  (2)(rl) 

3n  *  3 


Comparison  of  eqs.  ( 76 )  and  <  79 )  with  ( 38 )  shows  chat 


(1) 


(2), 


T0  ( a )  +  T_  ( s) 


T3(8) 


0 


a 


where  T  { s)  is  that  surfaoe  field  quantity  which  is  appropriate  when  B  and 
3  1 

B  are  electrically  connected.  If  T_(s)  has  already  been  computed,  it  is 
clearly  necessary  to  compute  only  one  of  Ty  '( s )  and  T  '  '( s ) . 

Let  us  now  return  to  the  potential  problem  set  forth  in  eqs.  ( 70 )  through 
( 73 ) .  As  regards  the  boundary  conditions  ( 70 )  and  ( 71 ) ,  an  exterior 
potential  satisfying  them  is 


(1) 


(2) 


2  *3 


(6 


where  V '  is  the  total  potential  considered  in  Section  2. 3.  Hence 


3$ 


3  n 


and  since  ( s )  is  given  as  the  solution  of  the  integral  equation  ( 35 ) ,  it  only 
.remains  to  specify  the  constants  7^  and  y^ . 

From  the  zero  charge  condition  ( 72 )  and  using  the  fact  that 
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// 


B1  °r  ®2 


dS  *  0  , 


we  have 


1  J  "  ‘3 
(1) 


pT^  1  s)ds  +  7  2  I  pT^(s)ds*-  J  pT9(s)  ds 


where  the  symbol  ( 1 )  below  the  integral  signs  shows  that  the  integrations  are 
carried  out  over  the  profile  of  the  portion  clone.  Similarly,  from  eq.  (73), 


J  pT3^(8)ds  +  y2  J  pTg  ^(s)ds 


?T2<  s)  ds 


where  the  integrations  are  over  the  profile  of  alone,  and  if  we  now  define 


CH  -  2»e  fpT3a)(8)d8  ,  C12*2^e  fpT^UJds  , 


C2l’2TC  f  ?hil){*)dB  *  C22  *  f  ^T3<2)(s)ds  , 


eqs.  ( 84 )  and  { 85 )  take  on  the  more  compact  form 


T1°U  +  Vl2  '  -2" 


p  T2(s  )  d  s  , 


(87  a) 
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71C21  ty2C22  *  '2t€ 


I 

(2) 


pT2(e)d8  . 


(37  b) 


It  will  be  observed  that  the  quantities  ,  etc.  all  have  the  dimensions  of 
capacity,  and  by  virtue  of  eqs.  ( 50 )  and  ( 81 ) , 


C11  +  C12  +  °21  +  C22  3  C 


(88) 


where  C  is  the  capacity  when  electrical  contact  is  maintained. 

Rather  than  solve  the  eqs.  ( 87 )  directly,  it  is  more  convenient  to  first 
eliminate  the  surface  Reid  quantity  T  {  s )  from  the  expressions.  That  this  is 
possible  can  be  shown  by  application  of  reciprocity  to  the  exterior  potential 
functions  z  -  V  ,  $  J  1  ^  and  $  ^  \  From  the  pair  z  -  V  and  $  J  1  ^ , 
we  have 


// 


$  o{1>  TTU-VjdS  « 


9n 


If 


B 


B 


Hence 


;/ 

B, 


If 


(1), 


T2(s)dS  *  JJ  zT2  '(s)dS  , 

B 


implying 


Similarly, 


/ 

(1) 

/ 


pT2(s)de  *  |  zpT3^(s)  ds 


l 

l 


(2), 


pT2(s)ds  *  I  zpTg  (s)ds 


(89) 


(90) 


(2)  *0 
and  wo  note  that  by  addition  of  the  last  two  equations  we  recover  eq.  ( 52 ) . 
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Finally,  from  the  function  pair 


§3(1)  and 


J  J  t}2\b)  dS  - 

B. 


implying 


P  T^2\b)  da 


T3(1)(s)dS 


B 


2 


pT3(1)(s)ds 


(91) 


i.e. 


(92) 


at*  expected. 

Using  eqe.  (89)  and  (90),  T  ( 8 )  can  be  eliminated  from  the  eqs.  (87) 

( 2 )  ^ 

and  if  we  also  eliminate  T  ( 8 )  using  eq.  (81 ) ,  we  obtain 

0 


(7l'VCll+V*t  ;  pT3(s)ds  *  -  2if  f  zc 

(i)  Jo 


zpT^  1  \s)  ds 


(7l“72)C21  +  72  2,rf  J  PT3(s)d8 

(2) 


J 


s 

l 


(  1  \ 

2  t  e  |  z  p  T3  (side 


-2»c 


s 


zpT3(s)ds  . 
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These  can  be  solved  to  give 


71  "  72 


“  f  f  zpTg(s)  de  J  pT3(e)dB-  f 
\*0  ( 1 )  Jq 


(1), 


pT3(s)ds-  |  zpT3'A,(s)dB  |  pT.^(s)ds 


s 

I, 


0 

(93) 


7  »  - 

2 


Jf  zpT3(s)ds  y  pT3U>(s)ds~  f  zpT3(1)(s)ds  f  pT3(1)(8)ds\ 

I/O  (1)  «*o  Jo  J 


where 


(94) 


s 

l 


/ 


(1), 


pT3(s)ds  I  pTg  (a)ds 


(1) 


s 

L 


pT3(1)(s)dB 


/ 


p  T  (  8  )  d  B  . 

*5 


(1) 


(95) 


We  can  now  proceed  to  the  calculation  of  P  .  If  we  write  this  quantity 

33 


as  F33  to  distinguish  it  from  the  Pg3  of  eq.  ( 55 )  for  B1  and  B2  in  electrical 
contact,  we  have,  from  eqs.  (8)  and  (83), 


33 


SI 


z  <T2(s)  +  71T3{1)(s)  +  72T3(2)(s)  >  dS 


B 


s  r 

w  jf  'p|t2(.)+  72T3(s)+(7l-T 


2>T^  1}(s)  >  ds  . 


( 9fi) 
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s 

-  2k  f 

n 


zp  (  Tg(  a )  +  7T,(s)  >  ds 


( see  eq,  46 )  where  7  is  given  by  eq.  ( 53 ) ,  and  thus 


P  »  I>  +  2ir  (y  -y)  J  zpT  (s)  ds  +  2 *  (yy  - t2>  I  zpTg  <s)ds 


33  33  "  '  ‘2  I  r  3 


Moreover,  from  eqs.  (53 ) ,  ( 94 )  and  ( 95 ) ,  after  some  manipulation. 


V1,  *  -(lfi*V~ 


pT3(1>(«)ds 


pT3(s)  ds 


which  enables  ue  to  write  eq.  ( 97 )  as 


p;,3  +  2.(7! -t2! 


( z  +  7)  pT  ^  ( 8 )  ds 
3 


The  factor  (7  -  7  )  is  defined  in  eq.  ( 93)  and  invoking  yet  again  the  expression 
(53 )  for  7  together  with  the  identity  (SI ),  we  have 


8 

' '  */ 


pTg(s)  ds  <7  I  pT,{s)ds+  I  zp’iy  (s)ds 
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2 ir  eA 


s 

/ 


(z  +  'V/PtJ  ^(s)  ds  , 


giving 


{Jt 


*33  •  p33  ■  TT^  r  U  +  T)pT„U><s>ds 


y 


(100) 


where  C  and  A  are  defined  in  eqs.  (50)  and  (95)  respectively. 


This  is  our  final  expression  for  P  .  Compared  to  the  situation  when 

JO 


B  and  B  are  electrically  connected,  the  only  additional  field  quantity  that 
1  2  /}) 

must  now  be  found  is  T  '  ( s ) ,  which  is  given  as  the  solution  of  the  integral 


equation  ( 76 ) ;  and  since  C/e  and  A  are  both  positive,  electrical  separation 
decreases  the  longitudinal  component  of  the  induced  electric  dipole  moment. 
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3.  AC0U8TICALLY  SOFT  OR  HARD  BODIES 

3.1  GENERAL  PROCEDURE 

Let  B  now  be  a  finite,  closed  acoustically  soft  or  hard  body  of  revolution 
about  the  z  axis  of  a  Cartesian  coordinate  system  (x,  y,  z).  It  It,  of  no 

concern  whether  3  is  disjoint  or  not.  A  plane  acoustic  wave  is  incident  and 

£ 

its  velocity  potential  is  written  as 


ikk.  r 


(101) 


A  g 

where  k  is  again  a  unit  vector  in  the  direction  of  propagation.  If  U  is  the 
scattered  field  that  is  produced,  then  U  8  satisfies 


(  V8  +  k2)  U8  »  0  In  V  . 


-|S-  -Iku8 

9  r 


0  as  r  -4  oo  , 


and  the  boundary  condition 


U®  «  -u1 


on  B 


if  B  is  soft,  or 


on  B 


(102) 


(103) 


(104) 


(105) 


if  B  is  hard.  Eqs.  ( 104)  and  ( 105 )  are  equivalent  to 


*  To  avoid  any  possible  confusion,  we  shall  henoeforth  refer  to  U  as  a  field. 


AhM1  MiD  nWl'tf 
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u  *  o 

on  B 

( lOfi) 

^  *  0 

on  B 

(107) 

8n 

respectively,  where  U  *  U  +  U  is  the  total  field. 

A  general  expression  for  U(r)  at  an  arbitrary  point  in  V  is  provided 
by  the  Helmholtz  representation: 


i  B 

where  R  «  jr-r’|  as  before.  For  sufficiently  small  k,  U  ,  U  and, 
hence,  U  can  be  expanded  as  power  series  in  ik  of  the  form 


Ui(  r ) 


(r)  > 


m»0 


( 109 ) 


and  when  these  are  inserted  into  eq.  ( 108 ) ,  the  coefficients  of  like  powers  of 
i  k  on  both  sides  of  the  equation  can  be  set  equal  to  give 


(£)+ irSteTT)!  //( 

B 


/  <  ,  \  -1-2  SR  ..  ,  ,  \ 

(m-f-l)R  > 


-  U^r’jjdS* 


(no) 


for  m  *  0,  1,  2, . . .  By  allowing  r  to  lie  on  B ,  an  integral  equation 
is  obtained  from  which  U  (i_)  can  be  found;  and  as  is  seen  by  substituting 
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the  power  series  for  U(  r )  into  eq.  ( 102 ) , 


o 

V"  u 


0 


ac 


V2U. 


3  0 


v2u 


m 


U 


m-z 


m  >  2  , 


(11 


showing  that  U.  ( r )  and  U  (r)  are  potential  functions,  but  U  (r)  is  not 

U  "  i  ^ 

unless  UQ(r)  =  0. 

In  the  far  zone  { r  —*  ao )  the  low  frequency  expansion  of  the  scattered 
field  deduced  from  eq.  ( 108 )  is 


US(r)~ 


ikr 
e _ 

4rr 


go  m 


1  1 

in*0  l»0 


(-1) 


m-f+ 1 


X 


.  r)  U^tr') 


+ 


( Kleinman,  1-965,  with  the  correction  of  a  sign  error),  provided  U  {  r ' ) 
is  taken  to  be  zero.  Our  objective  is  to  calculate  the  first  few  terme  in  this 
series. 


3.2  SQrV  BODIES 


We  now  specialise  the  above  results  to  the  case  of  a  soft  rotationally 
symmetric  body  illuminated  by  tho  plane  wave  ( 101 ),  and  seek  the  first  two 
terms  in  the  low  frequency  expansion  of  the  far  zone  scattered  field.  By 
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invoking  the  boundary  condition 


U  (r>  *  0 
m 


on  B, 


(113) 


*  0,  1,  2, _ _  we  have 


ikr  i'ffa 
U8(r)  U0^)dS'~ik 


0o(r')dSf  -ik  J  j[(r.  r')  ~ 

B 

U.  (rc)  !  dS*  +  °(k2)  >  » 

Bn1  1  “  J  J 


■£') 


(114) 


showing  that  only  the  potential  function*  UQ(r)  and  V^r)  are  required. 
From  eqs.  ( 101 )  and  ( 107 )  it  follows  that 


U*(r)  *  1  , 


Uj(r)  *  k.r 


(115) 


and  by  inserting  the  boundary  condition  { 113 )  into  { 110) ,  the  latter  becomes 


Um(r>  -  l£(r>  -  ^  Z(^T)l 

2*0 


// 


IT-  ut(i',dS’ 


which,  for  r  on  B,  reduces  to 


V1  (r) 
m 


m  mm 

j_  y  _i_  f  (k-1-1 

4a  «*  (m  '*)’•  J  J 

a  wf  W 


r~U.(r’)dS«  . 

dn!  *  ~ 


(116) 
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When  m*0#  eq.  (116)  gives 


..u-JJi 


rUn(r')dS‘  . 


9n  0 


This  is  identical  to  the  integral  equation  satisfied  by  the  potential  ( 
Section  2. 3 ,  and  henoe 


Vjr)  -  V  ( r)  , 


tT  u0<£)  -  VS> 


We  note  that 


s 

Jj?ruo(£>dS  m2’J 

B  ° 


)dS  *  2»  I  pT_(s)ds  »  — 

>3  € 


( see  eq.  50 ) ,  where  C  is  the  electrostatic  capacity. 
Front' eq.  (118)  with  m*l. 


A 

k.  r  . 
—  4 


B 


and  using  eq.  (  US),  the  left  hand  side  can  be  written  as 

,a  a*  .  ,a  <v  .  a  a  C 

( k  .  x )  p  woe  p  +  (k.y)psinp  +  (k,  z)z  -  ~ - 

a  5T  s 

Since  the  surfaoe  of  tho  body  Is  independent  of  jjt,  it  follows  that  U^( 
have  the  form 


imam 


(117) 


r)  of 


.  (118) 


(119) 


(120) 


r)  must 
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^(r)*  ^(k.x)co#0  +  (lc.y)ein0^j  U^1^(r)  +  (k.z)  U^'^tr)  0 


(121) 


(j> 


where  the  Individual  (jr),  j  *  1,  2  3,  satisfy 


■tjf 


s  ■ 


(122) 


B 


4r  J  J  R 


U^tr*)  dS-  . 


(123) 


B 


1  ■  It  Ilk  "t1  u‘3><£'»ds'  • 


(124) 


B 


Comparison  of  eqs.  { 124 )  and  ( 117  )  shows 


u/3)(r)  -  UQ(r)  -  V3(r)  .  *  T3<-J 


(125) 


(1) 


Similarly,  '  (r)  cos  0  is  identical  to  the  potential  V  (r  )  of  Section  2.2, 
implying 


(1) 


ly  (x>  "  V^p,  z) 


(126) 


so  that 


U,(1)(r)  »  Tj  ( s)  , 


an  1 


(127) 
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(2) 

and  U  (r)  is  identical  to  the  potential  V  (r)  of  Section  2.3,  so  that 


-k  ■  V>  • 


(12B) 


It  is  now  a  trivial  matter  to  evaluate  the  right  hand  side  of  eq.  ( 114). 
The  first  integral  is  clearly  C/c ,  and  the  second  can  be  written  as 

£*)  T3  i  ®  * )  “  ^jtfc.jbcoajJ’  +  (k.y)sin0^  T^s' ) 


+  (k.5)  1'2  (  8 f  )  -  T3(b*  «  >  dS’ 


]} 


«  2*  |  <  z)z'+ 


J  T:(s‘)  -  (k. 


z )  T2 ( 6 1 )  )  p ' d  s  1 


But 


s 

l 


s 

l 


P*T2( s»)ds’  »  I  z,p*T3<8»)  ds’ 


(52) 


and  hence  the  second  integral  on  the  right  hand  side  of  ( 114 )  is 


2  r 


8 

l  ( 


c 

4a  € 


y(r-£).  p'Tgts’lds' 


.  c  r_c_ 


y  ( r 


n  A 

-k).  zj 


where  >  is  as  defined  in  eqs.  ( 53 )  and  ( 54 ) , 
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The  low  frequency  expansion  of  the  far  field  is  therefore 


i  k  r  I  r 

+  °(k2>]  • 


(129) 


showing  that  a  knowledge  of  C  and  y  alone  is  sufficient  to  specify  the 
first  two  terms.  As  demonstrated  by  Van  Bladel  ( 1968 ) ,  a  similar  result 
obtains  even  for  a  body  which  is  not  rotationally  symmetric. 

3.3  HARD  BODIES 

The  final  case  to  be  considered  is  that  in  which  B  is  a  hard  rotationally 

symmetric  body.  The  boundary  condition  on  U  ( r  ) ,  m  *  0,  1,  2,  ....  , 

m- " 

is  then 


~  U  (r)  «  0 
3n  m  “ 


on  B 


(130) 


and  when  this  is  inserted  into  eq.  ( 112 ) ,  the  low  frequency  expansion  of  the 
far  zone  scattered  field  becomes 

U8(r)-^  (-1)  fik £  jin' .  r)U0(r’)dS’  +  k2f  Jj^r.r')  UQ(r' 

Lb  B 


-U^r'jJ  <n’. 


r)  dS *  -  ik3 


//jf(?.£'>2  »0CX*» 


~(r.  £* )  U^(jr_*)  +  U2(r')  (n'.rJdS*  +  Q(k*) 


( 131) 
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As  we  shall  see  later,  the  f'rst  term  0(k)  is  identically  zero,  and  we 
therefore  need  L’(  r ),  U.(r)  and  U,  (r )  to  compute  two  non-zero  terms 

U  —  1  U 

in  the  expansion. 

From  eq.  ( 110)  and  the  boundary  condition  ( 130) ,  an  expression  for 

U  ( r )  at  an  arbitrary  point  r  in  V  is 
m 


+  ft  2  (b7)!  ft  |)  IS1 

i«0  B 

(132) 


and  in  particular,  when  xn  *  0 , 


u0,£)  *  uo  <*>  +  ft 


// 


V£'>  ftr  ({)•*•  ■ 


(133) 


8  1 

Clearly  UQ  (r] )  *  UQ( r )  -  (jr)  is  an  exterior  potential  function 

r\  a  J 

and  —r —  U.(r)*  0  on  B  since  “7 — Un  (r)  *  0  .  In  addition, 

da  u  —  da  u  — 

s  •l 

U  (r)  vanishes  more  rapidly  than  r  as  r-*oo  since  there  is  no  term 

0(  k°)  present  in  the  expansion  ( 131),  and  hence 


U"  (r)  =  .  0 


(  134) 


implying 


UQ(r)  *  Uq(t)  -  1 


( 135 ) 


F  rom  eq.  ( 132 )  with  m  ■  1 ,  we  have 


U  { r '  »  U  ( r )  +  ~ 
1  "  1  -  4  r 
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which  can  lie  converted  into  an  integral  equation  for  r )  by  allowing 
r  to  approach  B .  Because  of  the  non-integrable  singularity  of  the  kernel 
for  r  on  B ,  it  is  necessary  to  apply  a  limiting  process,  and  if  a  bar 
across  an  integral  sign  is  again  used  to  denote  the  Cauchy  principal  value, 
we  obtain 


Uj  ( r)  * 


(r 


Ujtr* ) 


B 


(137) 


for  r  on  B,  where  uj(r)  is  given  by  eq.  (115).  In  terms  of  the 
cylindrical  polar  coordinates  (p ,  z), 


u|(r) 


.  x)  cos  0  +  (k.  y)  sinjj 


-} 


P 


a 

+  (k. z) z 


ami  since  the  surface  of  the  body  is  independent  of^,  it  follows  that  U1  ( r ) 
can  be  split  up  into  three  parte  each  of  which  has  the  0  dependence  of  that 
part  of  u*  ( r )  giving  rise  to  it.  In  particular,  on  die  arface. 


U^r) 


V  (s)  +  (k.  z )  V_ ( s ) 

4  o 

(138) 


where  V. ( a )  and  V_  ( s )  are  the  potentials  introduced  in  Sections  2. 4  and 
4  5 

2. 5  respectively  and  satisfying  the  integral  equations  ( 63 )  and  ( 67 ) . 

For  die  remaining  function  U  ( r)  an  expression  at  an  arbitrary  point 

M 

r  in  V  is  given  by  eq.  ( 132 )  with  n  ■  2  and  is 


Mr) 


(139) 
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where  (see  eqs.  101  and  109) 

U2  ^  *  2  *  (140) 

An  integral  equation  for  U2(,r )  can  be  obtained  by  allowing  r_  to  approach 

B ,  but  it  proves  unnecessary  to  determine  U  ( r )  explicitly  if  the  only  pur- 

3  * 

pose  is  to  calculate  the  term  0(  k  )  in  eq.  ( 131 ) . 

To  see  this  we  first  note  that  sinoe  UQS(r)  =  0  ,  the  eqs.  ( 111 )  imply 

y2  U2  *  °  ' 

8 

showing  that  U,  is  a  potential  function.  Moreover,  from  eq.  ( 139 ), 
s  1 

U2  *  u2  “  U2  iB  2111  exterior  Potential*  being  of  double-layer  type,  and 
since 


0  i  AAA 

~  U2  (r)  *  (k.  r)  (k.  n)  , 


(141) 


the  boundary  condition  on  U  (r)  is 

£ 


-k  U2(£>  - 


A  A  A 

- (k.  r) (k. n) 


(142) 


for  r  on  B.  clearly  depends  on  the  direction  of  incidence  as  well  as 
that  of  the  normal  to  the  surface,  and  in  principle  nine  separate  but  elementary 
potential  problems  must  be  solved  to  find  U®  .  In  tonne  of  these  potentials. 


**>  •  2  2 

i«i  j«i 

where,  for  convenience,  we  have  put 


kikjGij(^ 


(143) 
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i  I 
!£■  £ 


which  reduces  to 


J ^.(r'Hn’.xp-dS*  * J J#F|(r?)(n‘.xi)(r*.x.)dS' 


(148) 


B 


B 


when  the  boundary  conditions  ( 144 )  and  ( 146 )  are  employed.  Hence 

//u28(r')(n'.x4)ds,»2  Z  kikj  // 


G  (  r‘)(n' .  X|)dS' 


i  j 


B 


2  2 

i  j 


// 


kikj  JJV 

B 


*  (r'.xJdS' 


// 


A  A.  . 


F|(n' )  (k.  n* )  ( k.  r  * )  dS‘  , 


B 


implying 


JJ 


(r')(n!  r)dS* 


2//-, 


{ r  ')(r.  x^)(k.  n')(k.  r')  dS1 


B 


1  B 


and 


JJ 


U.  (r*)  (a* .  r)  dS* 


Z  jj  Y^r'Kr.  x£)  (k.  n‘)  (£.  r’ )  dS* 


B 


I  B 


+  2  //  rJdS* 


B 


(149) 
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This  integral  is  the  only  form  in  which  U  (x)  enters  the  far  field  expansion 
3  ^ 

through  terms  O  ( k  ) ,  and  sinoe  the  F^(  r ! )  are  known  by  virtue  of  the 
eqs.  ( 147 ) ,  the  integral  -an  be  computed  without  the  explicit  determination 
of  U2<x)  itself. 

We  are  now  in  a  position  to  evaluate  the  individual  terms  shown  in  eq.  ( 131 ) , 
Since  UQ(r)  »  1,  we  have 


// 


(n*.  r)  Un(r' >  dSf  ■  0  , 


(150) 


verifying  that  the  leading  term  in  the  far  field  expansion  is  ^  ),  and 


// 


(r .  r')  U  (r* )(n 


'.r)dS'«  *'jf 


*  \  A  O  * 


J V*{r.  r'Jdv' 


(151) 


where  VQ  is  the  volume  of  the  bo«ly.  Also,  from  eq.  (138), 


If 


U^r'Hn*.  r)dS*  »  r. 


2  v  e 

a 


.  A.*  .  A 

ere  or 1  cos  p '  x  +  cos  a  *  sin  0 '  y  -  sino ' 


x{  (k.  x)cos0'  +  (k.  y)sin  0*j  V^s* )  +  (k.  z)V5(  s' )^p* d  0*ds* 

A  /  A,  A  A.  A, A  a.  I  f 

*  r.  \  xl».x)  +  y(k.  y)  )x  j  p ! s’ )  cos  a'  ds ’ 
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where  a  *  is  the  angle  defined  in  eq.  (  59 ) .  Henco,  from  eqs.  ( 64 )  and  ( 68 ) , 


// 


U.  (r* )  (nf .  r)  dS’  *  k 


B 


■{< 


rM11  -  z(k.  z)(M 


11-M33>} 


(152) 


where  M1  ^  and  Mg3  are  the  elements  cf  the  magnetic  polarisability  tensor 
discussed  in  Sections  2. 4  and  2, 5  respectively.  As  we  have  previously  noted, 
for  a  body  of  revolution  Mg3  is  related  to  ( see  eq.  13). 

When  the  results  of  eqs.  ( 150)  through  ( 152)  are  substituted  into  eq. 

( 131 ),  the  low  frequency  expansion  of  the  far  zone  scattered  field  is  found  to  be 

US(r)~~  (k2  (?Mu-5(k.5)(Mu-M33)J  -V^j+  0(k3)j 

(153) 


where  the  actual  term  involving  k  is 


ikv 


JJ[i 

B 


( r .  r*  )^  U^(  r1)  -  ( r .  r' )  U  ^  r' )  +  U2(  r* )  (n!  r)  dS1 

— 


(154) 


Unfortunately,  the  evaluation  of  this  is  rather  a  messy  task. 


Since  UQ(r)  *  l  (see  eq.  135), 


// 


1  .A  .2  .  /A  I  A.  |  1  A 

-  (r.  r1)  Un(jr' )(n.  r)  aS  ■  -  r. 


J * f  n*  (r.  r* 


)2  dS* 


B 


B 
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//  (?.  r«)(£'.£)<k.  r')dS'  -  fr+(k.r)0  ■Iff  r 1  d  t1 


as  may  be  shown  by  analyses  similar  to  that  performed  above.  Hence 


j ' £  ( k.  r*)  (n*.  r)  -  (r .  £')  (n* .  k A  (k .  r ' 


)  dS* 


/// 


r’dr' 


( 159 ) 


which  cancels  the  contribution  { 155 )  of  the  first  term  in  the  integrand  of  ( 154 ) , 
The  complete  integral  (  154)  is  therefore 


ik3  j  F|(r')  jVr.x^  )(&.£')(  n'.k)  -  (k.Xj )(  r.  r ’)(n'.  r}^j  dS' 


( 161' ) 


I  B 


and  to  simplify  this  we  now  invoke  the  rotational  symmetry  of  the  body. 
From  eqs.  ( 147 )  and  ( 156 )  we  have 


( r * )  *  V^(  s' )  cos  01  , 


F  (  r 1 )  »  V  (s’ )  sin  0*  ,  F0(r’) 


V.(s')  . 
o 

(  161) 


When  these  are  inserted  into  ( 160)  and  the  azimuthal  integration  performed,, 
the  contribution  of  the  first  term  in  the  integrand  ir- 


jr  ( Ic.  z)' 


A  A  A 

(k. z) (r 


•■■oi 


p  ( z  au8  a  -  p  sin  a )  V  ( s )  d  s 


49 


013630-9-T 


+  T 


(r.  z)  j"l  -  {k. *)2J  J 


p  Vg{  8)  cos  o  da 


A  A  A  A  O 

2  jr  (r .  z)  (k.  z) 


s 

/ 


pzV  (s)  sine  ds. 

0 


The  contribution  of  die  second  term  in  the  integrand  of  ( 160)  differs  only  in 
having  r  and  k  interchanged,  and  when  the  two  are  subtracted,  the  final 

3 

expression  for  the  term  in  k  in  the  far  field  expression  ( 153 )  1b 


..3  .A  A  A 
•  k.  r  (k-r).  z 


(a  a  ,a  a,  .  a  C 

k.  r -(k.  z) fr.  z) >  I 


z)(r.z)^  |  p  (z  cos  or- p  sinn)  V4(s)  ds 


-(• 


,  A  A  .  A  A  . 

l  +  (k.  z)(r.  z) 


)/ 


p‘  V^(s)  cos  a-  ds 
o 


A  A,  A  A 

2  (k. z ) ( r , z ) 


8 

J. 


pz  V_(s)  sin  a  ds 
5 


{ 162 ) 


Although  this  is  only  the  second  non-zero  term  in  the  low  frequency  expansion, 

it  is  much  more  complicated  than  the  second  term  in  the  expansion  for  a  soft 

body.  The  surfaoe  field  quantities  involved  are  the  same  as  those  associated 

with  M  and  M  ,  but  there  is  now  no  simple  relationship  analogous  to  ( 13 ) 

11 

which  enables  us  to  dispense  with  v  ( 3 ) .  If  die  direction  of  incidence  or 

5 
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A  a  A  A 

observation  is  parallel  to  the  axis  of  symmetry,  i.e.  k  »  +  z  or  r  *  +  7.  , 
the  integral  containing  V^(  s )  disappears,  but  there  is  no  comparable  situation 
where  the  integrals  containing  V,.(  s)  are  absent  except  for  «he  special  case 
of  forward  scatter,  r  *  £  ,  when  the  entire  expression  ( 162 )  vanishes. 
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4.  THE  COMPUTATIONAL  TASK 

When  this  study  was  first  undertaken  the  main  objective  was  to  develop 
an  effective  program  for  computing  the  quantities  P  ,  P  and  M 

1 X  oo  X  1 

specifying  the  low  frequency  scattering  behavior  of  perfectly  conducting  rotation- 
ally  symmetric  bodies.  The  realisation  that  the  calculation  of  P  produces  as 

«3«3 

a  by  product  the  electrostatic  capacity  led  us  to  add  this  to  the  list  of  quantities 
considered,  but  it  was  only  later  that  the  question  of  acoustic  scattering  came 
up.  Since  the  first  two  terms  in  the  low  frequency  expansion  for  a  soft  body 
are  expressible  in  terms  of  C/e  and  y  ,  and  y  is  implicit  in  the  P 

33 

computation,  it  was  only  natural  to  add  this  to  our  list,  and  for  a  hart,  body 

the  first  term  involves  no  additional  work.  But  the  second  tern.,  ( 162 ),  is 

another  matter.  Ii>  particular,  it  requires  the  explicit  calculation  of  the 

surface  field  V  ( s )  that  had  hitherto  been  avoided  by  virtue  of  the  relation 
5  3 

{ 13 ) ,  and  even  if  this  were  done,  the  nature  of  the  k  term  is  almost  such 

ae  to  preclude  any  physical  understanding  of  the  data.  Foi  these  reasons  it  was 

decided  not  to  implement  the  computation  of  V  (  s)  and,  hence,  to  ignore  the 

second  term  ( 162 )  in  the  hard  body  expansion.  The  quantities  which  we  are  now 

left  with  are  all  ones  which  are  needed  for  the  electromagnetic  problem. 

4. 1.  INTEGRAL  EQUATIONS 

It  is  convenient  to  begin  by  listing  the  integral  equations  which  have  to 
be  solved  and  the  quantities  to  be  computed  from  their  solutions. 

Assuming  that  the  profile  p  *  p  (  z )  of  the  finite,  closed,  rotationally 
symmetric  body  has  been  specified  in  some  manner  and  its  volume  VQ  computed 
as  a  preliminary  step,  then: 
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{ i )  solve 


p*  K.T,  { a 1 >  ds 1 
x  1 


2  srp 


( 163  ) 


where  the  kernel  K  is  defined  in  eq.  ( 29 );  commute 

X 


f  f  P2T  (s)ds 

0  Jo 


( ii )  solve 


p»K0T2(8'>  ds*  *  2»z  , 


p  ’  KQ  I3  ( 8* )  d  b  ’  «  2 r  , 


(164) 


(165) 


(ir.fi) 


where  the  kernel  is  defined  in  eq.  ( 34 );  retain  the  option  to  print 
out  T^(  a);  compute 


C 


€ 


7  * 


pTj  (s)ds  , 


zp T^(  s  )ds  , 


(167) 


(168) 
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^21  , 

vo 


s 


zpT  (s)ds  -~~h~ 

0 


(169) 


( iii )  if  and  only  if  B  consists  of  two  separate  closed  parts  B  and  B  2 
solve 


s 

L 


2  x  r  on  B 


p*  K0T3(1)(s')ds’ 


;  (170) 


r  on  B. 


compute 


6  P 


33 


Jf 

volJo 


(1) 


(z+  7 )  P  T_  (s) ds 


(171) 


J 


(1) 


P  Tg  (  s)  ds  2 


(1) 


il 


p  T  (s)  ds 


where  the  symbol  ( 1 )  below  the  integral  sign  means  that  the  integration 


is  carried  out  over  the  profile  of  B  1  alone 


( iv )  solve 


1 {, 


(  s  * )  ^  p  cob  a'  0  2  + 


z )  sin  a  ’  -  p  '  cos  a  * 


«  x  ( 8  >  -  2p*| 


(172) 
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where  fl  ^ ,  Cl  and  a *  are  defined  in  eqs.  ( 61 ) ,  ( 62 )  and  ( 59 ) 
respectively  and  the  bar  across  the  integral  sign  denotes  the  Cauchy 
principal  value ;  compute 


8 


pV^(s)  cos  a  ds  . 


(173) 


We  therefore  have  four  ( five )  integral  equations  to  be  solved,  three 
(four)  being  of  the  first  kind  and  one  of  the  second,  and  five  ( six)  derived 
quantities  to  be  computed  from  their  solutions :  the  numbers  in  parentheses 
refer  to  the  unusual  situation  where  B  is  disjoint.  Before  attempting  this 
task,  there  are  certain  features  of  the  equations  to  be  examined. 


4.2  THE  KERNELS  AND  THEIR  SINGULARITIES 


The  kernels  and  of  the  integral  aquations  ( 163 ),  (  165),  (  166), 
( 170)  can  be  expressed  in  terms  of  complete  elliptic  integrals  of  the  first  and 
second  kinds. 

From  the  definition  of  R  given  in  eq.  ( 28 ) ,  we  have 


R 


r  2  2  'i  V2  2 

^(p  +  p*)  +  (z-z*)  J  ( 1  ~  m  sin  0) 


i/2 


where 


m 


lP-g-' 


(174) 


(175) 


(p  +  p  ’  )^  +  (z-z’)^ 


and 


i  <•'-*> 
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Hence 


\  1/2 


( 1  -  *n  sin  6 ) 


-1/2 


(177) 


and  when  this  is  substituted  into  the  definition  (  34 )  for  KQ ,  we  immediately 
obtain 


\  i/2 

m  i  ' 


K(m) 


(178) 


where 


t 

K{  m )  *  I  (1-m  sin2  8 )  i^2  d  8 

Jo 


(179) 


is  the  complete  elliptic  integral  of  the  first  kind  ( see,  for  example,  Abram owitz 
and  Stegun,  1964,  p.  590 ) . 

By  a  trivial  manipulation,  we  also  have 


COS  Ip 


£-1  i 

m  R 


R 

2  P  P  * 


(180) 


implying 


COS  Ip 


(mpp’)  1//2  (1-  m  sin“ 


2-1/2  9  i  fr> 

“0)  1/6 -{1-m  ahfer" 


(181) 


and  hence,  from  the  definition  { 29 )  of  K  , 


;1*(^7V/2  K(m)  -  E(n,)j  (.82) 


where 


E(m)  *  f  ( 1  -  m  sin2  8)  ^2  d0 


(183) 


56 


013(530-9-T 


is  the  complete  elliptic  integral  of  the  second  kind  ( ioc.  cit.  ). 

The  above  representations  of  KQ  and  K ,  are  exact.  Since  p  ,  z ,  p ' ,  z ' 
are  all  real  with  p ,  p  *  >  0 ,  it  can  be  verified  that  0  <  m  <  1 .  Over  this 
range  E(m)  is  a  finite  slowly-varying  function,  having  the  values  tt / 2  for 
m  *  0  and  unity  for  m  *  1 .  A  finite  polynomial  approximation  sufficient  for 

g 

computing  E(m)  with  an  error  of  less  than  2  x  10  is  given  in  Section  17.  3.30 
of  the  above  reference.  Through  the  first  three  terms  the  precise  expansion  is 
( Jahidr  nd  Emde,  1945 ) : 


E(m )  «  1  -  4mi  +  2mi^'  + 


(  184) 


m^  3  1  -  m  , 


„  (p-p')2  +  (z-g')2 

(p  +  p  '  )2  +  (  z-  z  '  )“ 


(185) 


( 180 ) 


r-  k  in 


*  M 

■  ~  in  — 

2 


(  187 ) 


We  observe  that  m^  *  0  if  and  only  if  p  r  *  p  ,  7.’*r. ,  that  is,  when  the  integration 
and  observation  points  coincide.  Fc”  an  integration  point  in  the  immediate 
vicinity  of  the  observation  point, 


(188) 


where  s  is  to  a  first  order  the  arc  length  between  the  points. 

The  elliptic  integral  K(m)  also  has  the  value  jt/2  for  m  *  0  but  becomes 
logarithmically  infinite  as  m  — ►!  .  A  finite  polynomial  approximation  sufficient 
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-8 

to  compute  K(  m )  with  an  error  of  less  than  2  x  10  is  given  in  Section 
17. 3. 34  of  Abramowitz  and  Stegun  ( 1984) ,  and  a  precise  expansion  through 
the  first  three  terms  is  ( Jafcnke  andEmde,  1945): 


K(m)  ■  P  +  2  “'ll-1  ~  4  mi  +  °^m12  *  V  )  • 


( 189) 


Because  of  the  infinity  of  K(m)  as  m  -*  1  (m^-4-  0),  KQ  and 
are  also  Infinite  in  this  limit,  but  their  behavior  in  the  vicinity  of  the  singularity 
is  easy  to  determine.  Using  ( 184)  and  ( 189 )  we  have 


(190) 


(191) 


showing  that  the  singularity  at  p  1  *  p ,  z  • »  z  is  an  integrable  one  in  each 

case.  The  contributions  of  the  singular  ( or  self  )  cells  to  the  integrals  in 

eqs.  (163),  (165),  (166),  and  ( 170)  are  therefore  finite  and  can  be 

analytically  approximated  as  follows. 

Consider  for  example  the  integral  equation  ( 165 ) .  If  the  self  ceil  in 

the  sampling  procedure  is  centered  on  s  *  s  (where  p  » p  )  and  is  of  arc 

□  □ 

length  As,  then 


p'K0T2(8»)ds* 


ln+2  ** 


sQ  -  -  As 


P’  K0T2(s,)ds‘ 
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in  +  7  A8 


-  P„T-<8n) 

n  E 


/u  ** 

, 


V18' 


s  -  —  As 
n  2 


~  T2  ^Sn) 


2*  AS 


"I  AS 


Bp 

4“  w  dsf 


and  hence 


I  ;‘K0T2(8')d8'  «  T2(8q) 


£  s  . 


It  is  desirable  to  retain  the  first  correction,  unity,  to  the  logarithmic  term  to 
ensure  the  necessary  accuracy  ehen  the  sampling  is  relatively  coarse  and  /  or 
p  is  small.  For  the  integral  equations  ( 166 )  and  ( 170)  the  results  differ 


from  the  above  only  in  having  TjU^)  and  '  UJ  respectively  in  place 

of  T  (s  );  and  for  the  integral  equntion  ( 163 ) : 

2  n 


s  +-  As 
_n  2 


p’K  T,  (s')ds*  «  P„Ti(sn) 


K,  ds’ 
1 
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giving 


/ 

self 


P*  K 1  Tx  (s')ds’ 


16p 

_ _i 

As 


As  . 


(193) 


For  the  integral  equation  (172  )  the  computation  of  the  kernel  is  a  more 
complicated  task  due  partly  to  the  presence  of  the  functions  Cl  ^  and  Ci^  . 
However,  these  also  can  be  expressed  in  terms  of  complete  elliptic  integrals, 
and  the  resulting  method  of  computation  is  much  less  time  consuming  than  a 
direct  numerical  evaluation  of  the  integral  expressions  for  Cl  ^  and  Cl  . 

The  definition  oE  SI  is  given  in  eq.  (61 ),  and  using  eqs.  ( 177  )  and 
( 181 ) ,  the  integrand  can  be  written  as 


I  \3/2 

cos  &  m  1  j  m 
r3  4  m  ^  p  p'  I 


{ 


( 1  -  ~)  (  1  -  m  sin"  9 ) 


-3/2 


-  ( 1  -  m  sin  9 ) 


f") 


from  which  we  have 


/' 


( 1  -  m  sin  0) 


-3/2 


d  9  -  K(  m ) 


194) 


To  evaluate  tho  remaining  integral,  differentiate  the  expression  ( 179  )  for  K(  m ) 
with  respect  to  m  to  get 


K’(m) 


1 

2 


s.‘n2  6 


(1  -  .i  sin"  & ) 


3/2 


d0 
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and  when  the  expression  ( 195 )  for  the  integral  is  substituted  into  this,  we  have 


1  f  m 


2  2  pp« 

m  v  r 


r  / 

/  i 


2m(l  -  — )  K’  (m)  -  ( 1  -  J  K(m'  +  E(m)  >  . 


(197) 


The  finite  polynomial  approximations  to  E  ( m )  and  K  ( m )  were  mentioned 
earlier,  and  in  particular,  for  the  latter, 

4  4  } 

K(  m )  *  ( a  +  a.  m  + . . . . -i- a .  m.  )  +  (b-+b,m1  + - +  b.m,  )  fn  - 

Oil  41  Uii  4  1  m. 


+  O  ( m.  3 ,  m  1  f  n  ~~  ) 
1  1  m^ 


(r*R) 


where  values  for  the  coefficients  a^  and  b  ^ ,  i  *  0, . . . . ,  4  are  given  in  Section 
17.3.34  of  Abramowitz  and  Stegun  (1964).  Since  d/dm  *  -d/-.m^,  it 
follows  that 


K’(m)*<  — ^-+  (b  -  a  )  +  (b  -  2a  )  m  -»■  (b  -  3a  )  m  “  +  (b  -  4a  )  m3 
|ml  1  1  2  2  1  3  3  1  4  4  1 


-(b  +  2b„m  +  3b„m  2+  4b.m.3  )  In— —  +  0(m  m  4ln  >  , 
1  21  31  41  1  1  m^ 


(  199) 


which  can  be  used  to  compute  K 1  ( m ) .  We  note  the  pole-like  behavior  of 
K'(  m )  when  in^  *  1-m  *0,  and  this  is  reflected  in  the  non-integrable 
singularity  of  the  kernel  of  eq.  ( 172 )  at  p '  *  p ,  z  *  *  z . 


4.3  THE  BODY  AND  ITS  VOLUME 


One  of  the  many  factors  motivating  the  present  rtudy  was  the  need,  to 
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compute  the  low  frequency  scattering  behavior  of  missile-like  targets.  These 
are  generally  rctationally  symmetric  bodies  (or  can  be  approximated  as  such 
to  an  accuracy  which  is  adequate  at  low  frequencies ) ,  and  are  often  made  up 
of  several  distinct  parts,  e.  g.  a  core  mated  to  a  cylinder  which  is  terminated 
in  a  spherical  cap.  Although  the  complete  profile  of  such  a  body  is  certainly 
not  an  analytic  curve,  each  individual  segment  has  a  relatively  simple  equation 
whose  form  can  be  used  to  advantage  in  the  numerical  process. 

It  is  therefore  assumed  that  the  profile  is  a  finite  piecewise  smooth 
curve  composed  of  etraighi  line  and  circular  arc  segments.  For  definiteness, 
the  number  of  segments  is  limited  to  15  or  less.  At  the  end  points  of  the 
rofile  where  it  intersects  the  z  axis  of  rotation  of  the  body,  p  =  0  (of  course ) , 
and  the  na  ure  of  the  program  is  such  that  segments  which  are  perpendicular 
to  the  z  axis  can  be  handled,  as  can  a  'disjoint'  body  having  two  separate 
parts  provided  each  portion  of  che  complete  profile  terminates  on  the  axis. 

Every  segment  contributes  to  the  total  volume  V(.,  which  can  be  found 
by  adding  the  individual  contributions  6  V,. .  In  certain  case  ,  a  volume 
contribution  can  be  negative  and  subtract  from  the  volume  attributable  to  the 
other  segments.  Where  this  occurs,  it  must,  be  noted  as  part  of  the  input 
specification  for  the  segment  in  question. 

In  the  following  we  list  the  input  specifications  of  circular  arc  ( Type?  1 
and  2  )  and  linear  (Type  3  )  segments,  and  give  exorcosions  for  the  cor¬ 
responding  volume  contributions  ( assumed  positive ).  The  segments  must 
be  described  sequentially  starting  at  the  intersection  of  Uie  left  hand  segment 
with  the  axis,  and  the  ordered  sequence  of  segments  defines  the  orofile  of 
tne  body.  .",i  some  cases  it  may  be  desirable  t>*  regard  a  single  linear  or 
curved  portion  of  the  profile  as  two  or  more  segments  to  permit  a  non-uniform 


r>3 
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spacing  of  the  sampling  points  or»r  the  whole. 

Type  1  Segment  (Circular  Arc .  Concave  Down) 

Specification: 


V  2  2 


',2“pU2) 


'K  degrees),  0  <  6  <  180 
volume  senbS 


O’  po 


If  6  is  the  angle  subtended  by  the  arc  at  the  center  of  curvature,  then 
the  radius  a  is 


a 


1 

2sta  | 


+ 


)‘ 


{ 200 ) 


Since  we  permit  the  specification  of  re-entrant  circular  arc  segments  we  do 
not  require  .  In  order  to  obtain  correct  results  for  both  standard 

(  z2  >  z1 )  and  re-entrant  segments  define  the  quantity 


d  - 


( 201 ) 
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Eq.  ( 200 )  gives  the  radius  a  of  the  type  2  segment,  but  the  coordinates 
( z qi  P0)  of  center  of  curvature  are  now 


i  f  o'! 

Z0  “  2  \^zl  +  z2  +  d(prp2)c0t2j 


2  p1  +  p2  +  d(z2"21)cot2  ‘ 


(  204 ) 


The  incremental  volume  of  the  type  2  segment  is 


6Vq  * 


r  j  !  V  Z1  ’  (P0  +  “2  "  I  <  u22  +  ulu2  +  u!  ’} 


•pt  Vp2'fOl'Ul(l’fi‘0,+  dS  8 


(205) 


Note  that  only  simple  sign  changes  distinguish  ( 202 )  from  ( 204 )  and 
( 203 )  from  ( 205 ) .  Relationship*  that  hold  for  both  type  1  and  type  2  segments 
may  be  derived  by  using  a  constant  5  defined  as  follows: 
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-1  type  1  segment 
1  type  2  segment 


( 206 ) 


The  center  of  curvature  (z Q,  pQ)  and  the  incremental  volume  6  VQ  for  circular 
arc  segments  of  types  1  and  2  are  then 


*0  *  2  {z1+22  +  5d(prt,2>00tf) 

p0  '  2  (f,l+p2  +  5d(z2_zl)00tfj 

2  ’  Z1  ’  (p0  +  “2  "  3  ( “l  +  Ul“2  +  U22  ’j  ‘  5p0  \_U2 ( P2  '  P0  > 


(  207  ) 


-  u.(p1-p0)  + da  6 


(203) 


where,  as  before. 


U2  *  V*0 


ui  *  zrzo 
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5.  NUMERICAL  C  -  .JTION  OF  THE  INTEGRAL  EQUATIONS 


p  M 

The  numerical  procedures  involved  in  finding  33  /  VQ  ,  11  /  V 

p 

and  where  appropriate,  33  /  VQ  are  quite  similar  to  these  required  for 

P  P 

11  /V  ,  and  it  is  therefore  sufficient  to  give  full  details  only  for  11  /  V ^ 


5.1  11  /V  COMPUTATION 


The  primary  task  is  the  solution  of  the  integral  equation  ( 163  )  for  the 

function  T  <  s)  and  this  entails  the  determination  of  a  sequence  of  values 
(i)  1 

T  ,  i  *  1,  2, _ _  N,  approximating  T^(s)  at  the  sampling  points 

s  a  s.  on  the  profile  p  *  p{  z ) .  For  this  purpose  the  profile  is  divided  into 
N  cells  Cj  of  arc  length  As  ^  and  midpoints  s.  corresponding  to  the 
coordinates  ( p^ ,  ) .  Within  each  cell  we  also  define  the  points  and 

8^+  where 


s  =s .  ~  ar  As, 

t  **  i  0  i 


■V  '  si  -  °o  isi 


(211) 


>v'ith  the  restriction 


0  <  -0  <  2 


(i) 


By  assuming  that  T  {  s }  has  the  constant  value  T  over  the  i  th  cell. 
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the  integral  on  the  left  hand  side  of  ( 163 )  can  be  evaluated  as  a  linecr  com¬ 
bination  of  the  whose  coefficients  depend  on  the  position  (p,  z )  of  the 

field  point,  leading  to  a  linear  system  of  N  equations  in  N  unknowns,  viz. 


(1) 


/ 


p'K^  ds' 


+  T 


(2) 


/ 


p,K1ds'  +  ... 


+  T 


(N) 


/ 


p 1 K1 ds ' 


2  up. 


N 


i  -  1,  2, 


N. 


(212) 


Hence,  the  system  to  be  solve'4  is 


A  t  *  b 


(213) 


where  t  ^  is  a  column  vector  with  elements 

t  »  T  (i) 

Cli  1 


J  *  1,  2 . .  N, 


(214) 


A  is  a  square  matrix  with  elements 


ij 


/ 


p,K1d8!  ,  ir  j*l,  2, . N,  (215) 


and  b  io  a  row  vector  with  elements 


bj  *  2  x  p  , 


j-1,  2,  ....  ,  N. 


(216) 


Increasing  the  complexity  of  the  quadrature  technique  used  to  evaluate 
the  integrals  f  will  generally  improve  the  accuracy  but  will  almost 


J 

C 


j 


?0 
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certainly  increase  the  computational  cost.  What  is  therefore  desired  is  the 
least  expensive  procedure  capable  of  giving  the  required  accuracy.  The  two 
simplest  approaches  are  t.o  integrate  first  and  second  order  approximations 
to  give  { i  y  j ) : 


aij '  j)  Asj 


(217) 


.j*  W1(V  ^j-Kl(iJ~>4  Pj  +  K1(i'^)J  +  wo(£>o)PjKl(i'  j) 


(218) 


respectively,  where  the  subscripts  j  -  and  j  +  correspond  to  the  points 
s.  and  s.,  of  eqs.  (211),  and  K.(i,  j)  is  the  kernel  defined  in  eq.  (  2S  ) 

j  -  j  +  1 

and  evaluated  at  the  points  (  p  . ,  z . ) ,  (p.  ,  z. ).  By  requiring  a  n  <  —  ,  w< 


3  J 


ensure 


that  the  sampling  points  s.  and  s.  do  not  coincide  with  the  end 

J-  3  + 


poin  ..  of  the  cell  ,  and  thereby  avoid  any  difficulty  in  the  computations 
of  n  and  Q  (see  eqs.  196  a*.  1 197 ) .  When 

A  4 


°om  iJj  • 


(219) 


eq.  (218 )  reduces  to  the  three-point  Gaussian  formula  for  which 


W  »  ***  VJ  ^  ~ 

0  9*1  18 


(220) 


With  this  choice  of  and  ,  the  advantages  of  eq.  (218)  vis-a-vis 
eq.  (217)  were  now  determined  by  computing  Pil/Vg  for  a  sphere  using 
various  values  of  N.  Fig.  1  shows  percent  accuracy  and  C.  P.  U.  time  versus 
N  for  8ach  integration  scheme.  It  is  apparent  that  for  a  given  expenditure  of 
C.  P.  U.  time  the  Gaussian  three-point  technique  is  much  more  accurate  than 


percent 

error 


3.0 


3.0 


0  10  20  30  40 


N 

■o 

Fig.  1:  Percent  error  xsd  C.  P.  U.  tims  of  *il/V^  cal¬ 
culation  for  a  aphere :  T  denote  trapezoidal  rule  eociputa- 
tioc  and  G  denotes  three-peiat  Giuaaiaii. 
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the  trapezoidal  method,  though  the  accuracies  of  both  are  severely  degraded 

if  N  is  too  small  (N<,5).  Since  the  Gaussian  scheme  with  N  *  10  produces 

an  accuracy  of  better  than  99. 8  percent  for  a  sphere,  there  is  nc  point  in 

going  to  a  more  complicated  procedure,  and  the  computer  program  was 

therefore  written  using  three-point  Gaussian  quadrature  to  determine  the 

matrix  elements  a . .  . 

ij 

In  summary,  the  integral  equation..  ( 163 )  is  solved  by  conversion  to 
tbs  matrix  system  (213)  in  which 


*ij*  18  (vKl<1',->  +W‘,i+>j+J9  p)Kl(l*J)  Asj 


i,  j  *  1,  2, - -  N;  i  f  j 


(221) 


(l)  P 

Having  determined  the  sampled  values  =*  T  (s  ),  ll/VQ  is 

computed  from  eq.  ( 164 )  by  integration  over  each  segment  of  the  profile 
using  a  second  order  integration  procedure  (  subroutine  INTEG ,  described 
in  the  Appendix ) . 


5.2  P33/Vq  COMPUTATION 


Th  .  point  sampling  method  of  solution  of  the  integral  equations  ( 165) 
ani  ( 16o ,  requires  us  to  .find  the  sequences  T  ^  3  T  (s.)  and  T  ^  *  T  { u. ), 

b  £  1  #,\  o  *3  1 

i  *  i,  2, - -  N,  from  these  equations.  To  determine  the  T  ,  choose 

0^,  w  and  w,  in  accordance  with  (219)  and  (220)  and  thence  solve  the 
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matrix  system  A 1 2  *  b  where 


t  *  T 
c2i  2 


(i) 


i  *  1,  2, _ ,  N 


(222) 


bi  *  2rzi 


and 


ij 


[fa  {vKo(i'>"> 


+  pJ+Ko(i-i+)  r  ?  <>}*■ 


ii 


s 

In  | 

\ 

iA8i  /  J 

As 


i  1 


0(i.j>  I  AS 

i,  j  *  1,  2, . N;  i  yf  j  ,  (223) 

i  *  1,  2,  .....  N. 


(i) 


The  T  are  similarly  determined  by  solving  the  matrix  system  At  ~  b 

O  O 


where  the  elements  a  ^  are  again  given  by  ( 223 ) ,  but 


t  *  T 
1 3i  3 


(i) 


i  *  1,  2,  ....  ,  N  . 


(224) 


b  ■  2r 

i 


C  p 

The  quantities  /e,  y  and  33 /V^  defined  in  eqs.  (167),  (168)  and  (169) 


respectively  are  computed  using  the  same  second  order  integration  procedure 


employed  in  calculating  *11 /V  . 


If  the  body  profile  nonsisia  of  two  discrete  partst  it  is  also  necessary  to 
solve  die  integral  equation  ( 170).  The  corresponding  matrix  system  is 
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almost  identical  to  that  in  ( 224 ) ,  and  from  the  sampled  values  T  ^  ^  (  s  . ) 

O  i 

ru 

6  P  P 

and  T^Sj),  33 /VQ  171 )  is  computed  and,  hence,  33/V() 


5.3  11  /VQ  COMPUTATION 


The  basic  approach  is  similar  to  the  above  in  spite  of  the  more  complicated 

integral  equation  ( 172 )  that  must  now  be  solved.  The  matrix  equation  for  the 

sampled  values  V,( s.)  *  V  ^  is  Av  »  b  where 
4  i  4  4 


v  *  V 
4i  4 


b  *  2irp . 
i  l 


i  3  1  ')  K 

*  *  p  •*!  •  •  •  *  t 


(225) 


'  [fa 


)  +  p.  +  f(i,j  +  A  +  ~  ASj 


i,  j  =  1,  2, - -  N;  i  f  j  , 


( 220 ) 


s .  +  r  A  s 


"  ‘  f 


p  !  f  (8  ' )  ds  ' 


i  3  1,2, . N 


v  5  Asi 


in  which 


r  r 

f( i ,  j)  *  p  cos  q.  Q  (i,  j )  +  (?..  -  7.  )  eincr.  -  p.cos  o  . )  n  (i ,  j ) 

1  J  1  L  J  1  J  J  ]J  1  J 


i,  j  *  1,  2, _ _  N;  i  j  j  . 


( 227 ) 
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We  observe  that  the  computation  of  each  diagonal  element  of  A  requires 
the  numerical  evaluation  of  a  Cauchv  principal  value  (denoted  by  the  bar  across 
the  integral  sign  in  the  above  exp;  ession  for  a  ) .  As  an  approximation  to 

this  principal  value,  we  remove  from  the  cell  C  a  slice  defined  by  the  interval 

1  1  1 
( s  -  7*  /3  A  s  ,  8  +  —  /3  As  /  where  (3 ,  0  <  [3  <  1 ,  is  the  fractional 

1  JL  i  1  4_l  1 

exclusion;  (3*1  implies  no  exclusions,  i.e.  that  the  principal  value  is  not 


taken. 


We  now  have 


s  -  ^  A  s 
"i  2  i 


V  2  ASi 


p*  f(s')  ds'  - 


p’f(s')  ds' 


s  -“As 
i  2  i 


8i+  f  isi 


(228) 


and  these  integrals  are  also  computed  using  three-point  Gaussian  quadrature. 


Defining 


8i  “  4  *si 


8il  ”  Si2  ”  2  "o'1-131  A81 


13  *  812+2°0*1’*!*A::: 


Si5  *  8i  +  4  U+|S)  A8i 


(  229  ) 


8H  *  St5  -  2  °0U‘  A£i 


8 1 6  ”  s15  +  Ja0(I-P>  A8i 
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v  3  obtain 


5  ( 


,i  -  »  -  I(1-8i  |  18  <  pUf(i'1l)  +  pi3‘<1',3)  +  pi4,<i-i4) 

*N  "\  I 

+  p  ,  e  f  ( i ,  i  „ )  )  +  “  /p._f(i,i„)+p  f(i,i)  \  I  A  s  . 


(230: 


Equations  (225)  through  (230)  completely  describe  a  system  of  N  linear 


equations  in  N  uukncwns  V^,  i  =  i,  2, _ _  N.  Their  solution  and  sub¬ 
sequent  integration  of  the  V  .  according  to  eq.  ( 173 )  yield  M  ,  . 

4,  11/ V0 

Experiments  were  performed  to  find  an  appropriate  value  for  the 
fractional  exclusion  ft.  Ab  an  example,  for  a  sphere ^Mn  /\'0  *  1,5^ 
with  N  *  20,  the  data  in  Table  1  were  computed.  If  we  exclude  the  fortuitous  (?) 


Table  i 


p 

Mn/v0 

percent 

error 

1.0 

1.480 

-1.33 

C  1 

1.  516 

j..  07 

0.01 

1 . 503 

0. 18 

0.  001 

i.  501 

0.  08 

error  zero  occurring  for  ft  somewhere  in  the  range  (.*.  1  <  ft  <  i ,  these 
data  indicate  that  the  choice  ft  =»  u.  001  is  sufficient  to  keep  the  error  It* 
than  0. 1  percent. 
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Increasing  the  number  N  of  points  at  which  the  surface  is  sampled 

will  generally  increase  the  accuracy  of  computation,  but  since  the  number 

2 

of  matrix  elements  iicreases  as  N  and  the  cost  of  a  linear  system 

3 

solution  increases  roughly  as  N  ,  this  improvement  is  obtained  at  the 
expense  of  an  increase  in  computation  cost.  Unfortunately,  there  is  no 
rule  for  specifying  the  minimum  value  of  N  sufficient  for  a  given  accuracy, 
and  the  information  which  follows  is  based  only  on  our  experience  in  using 
the  program. 

P 

The  results  in  Fig.  1  and  Tablt  1  show  that  for  a  sphere  il/v  and 

M  ^  ^ 

11  /V  are  acouvataly  determined  with  N  as  amc.ll  as  20,  and  this  is  also 

jp 

true  oi  33  /  VQ  .  On  the  other  hand,  if  the  body  has  a  discontinuity  in 

dp/dz  iying  off  the  axis,  it  appears  necessary  to  increase  N  to  50  or 

P  i 

more  to  maintain  the  same  accuracy  (error  £  0.5  percent)  in  the  11  /V 

P  ^ 

and  33  /  computations.  This  i*  illustrated  by  the  results  in  Table  2  for 


Table  2 


N 

Ni 

*2 

pu/v 

11 

7 

4 

2.752 

20 

15 

5 

2.801 

40 

30 

10 

2.872 

70 

50 

20 

2.888 

Note  :  W  is  the  number  of  sampling  points  oa  the 
generator  of  the  coee  ( linear  segment)  and  is 
the  number  on  the  (half,1  base  { circular  arc  segment) . 
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a  rounded  cone  with  halt  angle  15  .  The  small  but  not  negligible  ( 0. 58 

p 

percent)  change  in  11  / ac  M  increases  from  40  to  ?0  suggests  that 
such  large  values  of  N  may  be  essential  for  bodies  such  as  tins  for  which 
T^(s),  T?(s)  and  T^(s)  have  infinities  at  one  or  more  points  on  the 
profile. 

M 

For  the  same  rounded  cone,  the  results  tor  11  /  V  are  given  in 
Table  3.  Since  an  increase  in  N  front  17  to  35  produces  only  an  insignifi- 


Table  3 


N 


N 


1 


N, 


M 


11 /V, 


17 

35 


10 

25 


7 

10 


1,880 
1  678 


M 

cant  change  in  11  /  V_  ,  the  choice  N  3  20  is  now  adequate.  Observe 

°  M 

that  the  surface  field  V  .  (  s  )  associated  with  11  V,.  does  not  become 

4  '  0 

infinite  at  a  discontinuity  in  dp  /  dz,  and  this  is  undoubtedly  the  reason 

why  in  many  cases  a  email  value  of  N  new  produces  the  same  accuracy  as 

P  p 

does  a  much  larger  value  in  the  11/ V„  and  33/  V  computations. 

0  0 

Nc  attempt  has  been  made  to  exploit  this  finding  in  the  general  program. 

When  treating  bodies  composed  of  several  segments,  a  strategy’  which 
has  proved  successful  Is  to  divide  all  segments  into  cells  of  approximately 
equal  length.  This  serves  to  fix  the  allocation  of  any  given  nu  tber  X  of 
sampling  points  among  the  various  segments.  Tests  sc  far  per  >rmed  have 
not  conclusively  shown  the  advantages  of  dividing  a  single  segment  into  two 
or  moire  smaller  segments  so  as  to  c-’feet  a  non— .inifomi  sampling.  It  is, 
however,  believed  thst  such  a  sub-division  may,  for  a  given  N ,  improve 
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6.  CONCLUDING  REMARKS 

We  have  here  considered  the  low  frequency  scattering  of  electiomagnetic 
and  acousf  ;  waves  by  axially  symmetric  bodies.  By  concentrating  on  certain 
quantities  such  as  the  normalised  components  of  the  induced  electric  and 
magnetic  dipole  moments,  we  have  shown  how  it  is  possible  to  arrive  at 
rather  elegant  expressions  for  the  far  zone  scattered  field  in  terms  oi 
quantities  which  are  functions  only  of  the  geometry  of  the  body  Each  such 
quantity  is  expressible  as  a  weighted  integral  of  an  elementary  potential 
function  which  can  be  found  by  solving  an  integral  equation. 

A  computer  program  has  been  written  to  solve  these  equations  by  the 
moment  method  and  to  calculate  the  dipole  moments,  the  electrostatic 
capacity,  and  a  further  parameter  y  related  to  the  capacity.  Any  body 
can  be  treated  whose  profile  is  made  up  of  straight  line  and  circular  arc 
segments  and  it  is  even  possible  to  have  two  distinct  bodies  with  or  without 
an  electrical  connection  between  them.  Although  no  Berious  attempt  has 
been  made  to  optimise  th<"  program,  only  a  few  seconds  are  required  to 
compute  ail  of  the  aDove  quantities  to  an  accuracy  of  better  than  one  half 
percent. 

We  have  already  used  the  program  to  compute  the  scattering  from  a 

variety  of  shapes,  and  it  may  be  helpful  to  list  some  ci  the  results  obtained 

so  far.  Data  for  a  rounded  cone  cons,  g  of  the  intersection  of  a  cone 

of  half  angle  0  with  a  sphere  centered  on  the  apex  are  given  in  Table  4. 

f/w  is  the  length -to -width  ratio  of  the  body.  For  0  v  90°,  the  values 
P  P 

of  11  /  Vq  and  33  /V^.  are  quite  similar  to  those  previously  com¬ 
puted  by  Senior  ( 1971 )  using  a  mode  matching  method,  but  since  ^  11  /  V 
showed  significant  discrepancies,  this  quantity  was  determined  for  a  variety 
of  0.  Detailed  checking  has  confirmed  that  tre  present  data  .-ire  accurate 
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Table  4 :  Rounded  Cones 


(deg.  ) 

l/w 

vo 

Pn/v0 

P33/V0 

M“/vo 

C/(e/fw 

3 

9.554 

0.  002870 

1.884 

7.5 

3.837 

0.01792 

1.813 

15 

1.932 

0.07137 

2.865 

c.  147 

1.678 

30 

1.000 

C.  2806 

3.664 

3. 4,94 

1.484 

5.  406 

45 

0.7071 

0.  6134 

1.366 

60 

0. 5774 

1.047 

4.520 

1.931 

1.312 

6.386 

80 

0. 5321 

1.731 

1.334 

90 

0. 5000 

2.094 

4.428 

2. 184 

1.373 

93 

0. 5262 

2.204 

4.368 

2.242 

1.386 

7.428 

99.2 

0. 5799 

2.429 

4.261 

2.372 

1.416 

7.303 

108 

0.6545 

2.742 

4,071 

2.553 

1.458 

7. 123 

120 

0.  7590 

3.1' 2 

3.789 

2.769 

1.507 

6.  889 

140 

0. 8830 

3.699 

3.370 

3.006 

1.  547 

f>.  586 

151.7 

0. 9402 

3,938 

3.187 

3.042 

1.540 

6.441 

180 

1.0 

4. 189 

3.0 

3.0 

i.  5 

6.  283 
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to  three  significant  figures.  Whan  9  *  90°  the  body  is  a  hi  li sphere 

P  P 

for  which  precise  values  of  11  /  V  and  33  /  V  are  available: 

P  p  u  u 

11  /V  •  4.430 _ _  33  /Vq  *  2.139....  ( Schiffer  and  Szego, 

1949,  p.  152) .  The  corresponding  values  in  Table  4  are  within  0.  2  per¬ 
cent  of  these.  For  6  >  90°  the  cone  is  a  re-entrant  one,  i.o.  a  sphere 

Q 

with  a  conical  region  removed,  and  when  9  *  180  the  bxly  is  a  sphere 
for  which  exact  data  are  also  known. 

Results  for  ogives  and  symmetrical  lenses  whose  arcs  subtend  an 
angle  9  at  their  centers  of  curvature  are  shown  in  Table  5  .  The  trans¬ 
itional  shape  is  a  sphere  for  which  9  *  180°. 

To  illustrate  the  computations  when  two  bodies  are  present.  Table 
6  gives  data  fGr  two  identical  spheres  separated  by  a  distance  i  d  where 
d  is  the  sphere  diameter.  When  the  two  c  oheres  are  touching  ( e  *  0) 

It  is  known  that  **33  /  *  8/3  ( Schiffer  and  Szego,  1949,  p.  154); 

the  ratio  deduced  from  Table  6  is  2.  678,  which  is  within  0.  4  percent  of 

P  M  p 

the  exact  value.  As  e  increases,  11  /V  .  11 /V.  and  33/ V,. 

u  0  u 

rapidly  approach  the  values  appropriate  to  a  single  sphere  in  isolation. 

p 

33  /Vq  ,  on  the  other  hand,  is  proportional  to  the  axial  cc  nponent  of 
the  induced  electric  dipole  moment  for  two  spheres  which  are  electrically 
connecteu  by  an  infinitesimal  wire,  and  with  increasing  e  this  increases 
indefinitely,  as  expected  (  Kleinman  and  Senior,  1972).  The  same  is  true 
of  C/.  x  .  The  parameter  7  has  also  been  included  in  Table  6,  and 

Vq/ 

since  its  exact  value  can  be  shown  to  be  -( 1  +  e/2  ) ,  the  accuracy  of 


computation  can  be  judged. 
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Table  5  :  Ogives  and  Lenses 


Shape 

0(deg.  ) 

if  w 

v„  P 

n/v0 

P33/V0 

o 

> 

r-t 

s» 

C /(c'Jt w  ) 

/ 

22.? 

10.  C? 

0.  004146 

2.089 

49.88 

1.943 

6.128 

36 

6.314 

0.  01053 

2.100 

24.15 

1.910 

5. 696 

ogi 

ve 

56 

4.011 

0. 02650 

2.189 

12.57 

1.843 

5. 451 

88 

2.475 

0. 07148 

2.363 

6.778 

1.739 

5.424 

132 

1.540 

0. 1966 

2.647 

4.136 

1.611 

.”>.  696 

' 

f 

150 

1,303 

0. 2847 

2.775 

3. 595 

1. 564 

5.880 

sphere 

180 

1 

0. 5236 

3 

o 

o 

1.5 

6.283 

j 

\ 

64.4 

0. 6297 

1.212 

3.779 

2.161 

1.368 

6.570 

43.6 

0. 4000 

2. 586 

5.182 

1.674 

1.252 

7.375 

lei 

18 

28 

0. 2493 

6.448 

7.649 

1.390 

8.758 

17.2 

0. 1512 

17.  30 

11.88 

1.225 

1.098 

10.  759 

t 

11.4 

0. 09981 

39.55 

17.49 

1.144 

1.  061 

12.982 
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Table  6 :  Two  Spheres 


e 

Pa/v0 

P33/vo 

P33/Vq 

Mu/v0 

c/(^>  7 

0 

2. 702 

7.237 

7.237 

1.605 

6.153 

-1.0004 

0.005 

2.706 

7.280 

5.021 

1.586 

6.151 

-1.0029 

0.01 

2.709 

7. 319 

4.800 

6. 150 

1 

►-* 

o 

o 

4* 

0.02 

2.715 

7.402 

4.557 

1.607 

6. 148 

-1.0100 

0.  05 

2.732 

7.655 

4.  ?.i0 

1.592 

6. 141 

-1. 0250 

0. 1 

2.759 

8.026 

3.922 

1.579 

6. 130 

-1.0503 

0.  5 

2.391 

12.02 

..299 

1.528 

6.  015 

-1.2502 

.1.0 

2.950 

18.19 

3. 142 

1.511 

5.822 

-1.5001 

5.0 

2.994 

120.6 

3.045 

1.500 

4.304 

-3.5001 

10.0 

3.  002 

383.1 

3.046 

1.501 

3.470 

-6. 0000 
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APPENDIX  :  THE  COMPUTER  PROGRAM 

P  1’ 

The  program  computes  11 /VQ  ,  C/£,  33/ V  , 

A# 

M  P 

11  /Vq  and,  where  appropriate,  33 /Vq  ,  and  consists  of  a  main 

program  and  six  subroutines. 

A.  1  DATA  SET 

A  data  set  is  made  up  of  one  control  card  and  a  Dumber  of  segment 
specification  cards,  one  for  each  segment  (or  sub-segment)  of  the  profile. 
The  segment  specifications  conform  to  the  convention  stated  in  Section  4. 3 . 

Control  Card 

Columns  Description 

1  The  number  ( 1  or  2 )  of  bodies. 

3-4  Two  digit  integer  ( right  justified ):  the 

number  of  segments  on  the  first  body 
( the  body  to  the  left ) .  When  there  is 
only  OLe  body,  use  these  columns, 

C-7  Same  as  columns  3-4,  but  for  body  to 

the  right. 

9  A  printing  key : 

p 

1 :  print  T^  from  33  /  com¬ 
putation. 

0  or  blank :  do  not  print  T  . 

o 
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s 


11 


13 


21-30 


A  computation  key  ( 0 ,  blank  or  1 ) 

o 

1 :  suppresses  commutation  of  11/  Vq  > 
C/E.  T,  P33/V0. 


A  computation  key  { 0 ,  blank  or  1 ) 


M  , 

1 :  suppresses  computation  of  11/ 


A  real  number:  the  fractional  exclusion  ft. 
If  these  columns  are  blank,  ft  defaults  to 
0.001 . 


Segment  Specification  Card 
Columns 

1-2 


4 


r 


11-20,  21-30 


Description 

Two  digit  integer  (  right  justified) :  the 
numoer  of  sampling  points  or  cells  on  the 
segment. 

Segment  type  key : 

1 :  circular  arc,  concave  down 
2  :  circular  arc,  concave  up 
3 :  linear . 

Volume  sense: 

+  or  blank  :  additive  volume 
-  :  subtractive  volume  . 

Two  real  numbers  :  respectively,  the  end 

coordinates  z  and  z  of  the  segment. 

1 
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31-40,  41-50  Two  real  numbers :  respectively,  the  end 

coordinates  p^  and  p„  of  the  segment. 

51-60  A  real  number:  for  circular  arcs,  the 

included  angle  in  degrees. 

There  are  the  following  restrictions  : 

( i)  the  total  number  of  segments  must  not  exceed  15, 
and  ( ii )  the  total  number  of  cells  over  all  segments  must  not  exceed  80. 

The  profile  is  specified  i.i  tbs  direction  of  increasing  profile-length,  beginning 
at  its  left-hand  intersection  with  the  z-axis  and  ending  at  its  right-hand  intersec¬ 
tion  with  the  z-axio.  Re-entrant  segments  are  permitted,  allowing  z^  >  z^  . 
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A.  2  MAIN  PROGRAM 


The  main  program  reads  and  prints  data  and  supervises  all  computations. 
A  rough  flow  chart  showing  the  interaction  of  the  subroutines  is  given  below 
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RF-M  API  1  !  HO,  HO)  ,AP33(RO,HO  1  ,AM1  1  {><(),  HO)  ,  X  <  MO)  trt(hU)  . 

4  7  F  P  (  ?  )  ,RHOFP<?  )  .ST3v4)  ,T3( RO),Mll  ,T1  (HO) 

INTFGFk  NUMPTSUS)  ♦  PI.OS  /  »  +  •  /  ,m.  /  •  »  /  ,  1  NQX  (  7  1 

COMMON  RHO(RO,9  )  ,7  (80,9)  »  ARC  <801 ,C(R0,9)  ,  S  (  HO ,  9  )  /  Slit.  /  I  *-»  S  C  80  i 
DATA  MIN.TWOPI ,PI ,WO,Wl / ,  4 . ?H 31 8 5 , 3 . 14 1 593 , ,444444*. ,?7  7  1111/ 
3?  WFAO! 5,34, FND=999  )NHOD»NSl  ,  MS?  ,  !  PR  I  NT  ,  KFYP 11  ,  *  F  YM 1 1  ,  PR 
34  FORMA  T(  I1,1X,?(I?,1X),3!I1  f  IX  ;  ,4X , FI ?.7 1 
WR !  T  F (  A,4)NH0r,,NSl 

a  FORMA  T(  '  1***  KFOIMNir-T,  OF  DATA  S  FT  :  ‘  /  •  O  •  ,  5X  ,  *  MOD  1 1:  S  '  ,  A  X  ,  •  -  •  ,  I  ?  / 

3  •  •  ,  5X  ,  •  SfcGMFNf  5s  :  *  /  '  •,5X,‘RnOY  h  \  *  ,  5X  ,*=',!  ?  ) 

lF'MPOO  .LF.  0  .OR.  "ROD  .GT.  ?)  GO  10  <190 
NSFGS  =  NS1+0'v? 

I  F  (  MR OO  .  FO.  ?)  WR  I  TF  <4  .1  INS? 

I  FORMAT!*  *  ,  5X  »  •  RODY  »?•  ,bX, '  =  '  ,  1?  ) 

1F(NR0I)  . FO.  1)  GO  lO  ]0 

IF  (MSI  .  I.  F .  O  .OR.  f«is?  .LF.  0)  GO  TO  <<40 
10  TF(NXPGS  .IF.  A  .OR.  NSFGS  .GJ.  IS)  GO  lO  990 
WR  1  TM  4,  35?  1  I  PR  INT  ,KFYPU  «KFYM1  1 

3S?  FORMAT!'  IX, 'PRINT  K  F  Y  '  ,  3X  ,  •  =  '  ,  I  ?  /  •  '.Sy.ifllMP  t^Y  P  •  •  5  •  •  ,  I  ?  / 

4  •  ',SX,'COMP  KFY  M'*S  =  ',!>') 

IF!  t  p  ,I.F.  o.  .OR.  FR  .(,1.  1.)  FR  =  .0.)1 
I  F  (  K  F  YM)  l  .Fo.  01  WRITF(4,5)FP 
S  FORMAT!  •  •  ,*>X, 'FXr.I.MSlOM  :  '  ,  F  /  ,/«  ) 

1F!RFYP11  .NF.  0  .AND.  KFYMIJ  .OF.  0)  (U)  TO  4'R. 
v.  =  0 
NCI  =  0 
V0=0.0 

oil  11  1  =  1,  NSFGS 

-FAD!  5,1?)  NOMPTS (  l  )  ,  I  T YP  ,  I  S I r,N-  Z FP ,RMOFP ,  ThF T  A 
1 ?  RORMAl (1 ?, 1 X, II  , 1 X,A1 ,4X,SF10.7  1 

IFJMIMPTSU)  .IF.  O  .OR.  IT'P  ,I.F.  O  .OR.  I  TYP  .Gi.  X)  )i  090 

IF!  I  S I  GN  .FO.  ML  I  IMGN=PI.iiS 

'-JR  I  TF  (  4,1  3)  I  ,0'iMPTS!  1  )  ,  ITYP,  lSIf,N,7Fp,RH!iFp 

13  FORMA  S{  'OShGMFNT  **•  ,!?,•:«/  *  '  ,  5X  »  '  C  F I  I  S', 7X, '  =  ',!?/• 

A  tTYpt-  KFY'  ,4X  ,  *=  •  ,  I  ?/ •  '  ,*>X,  •  VOUIMF  SFNSF=  «,A4/«  '.<•>. 

R  »7-C.0:iRDINATF  FNO  POINTS  =(•  ,  F  )?.  7  F  I  ?./,')*/ •  '.Sx, 

5  1  P HO— COORn  I N  ATF  FNl)  PO I  NT  S=  !  •  ,  F  1  ?  .  7  ,  '  ,  «  .  H  ">  .  !  ,  •  )  •  ) 

I  F  (  I  TYP  .NF.  3)  WR1TM4,1A>  THF  T  A 

14  f  ORMA 1 ( •  • ,5X»  ‘ TMFTA  !OFG)  = ' , F 1 0 . S ) 

IF(KFYM11  .FO.  1)  FR=1.0 

f-‘  ~  Vi 

IMF  TAaI’I*THFl  A/IHO. 

MrM^-MHMpT  S  (  l  ) 

l  F  (  I  .IF.  NS))  Nf  1  =NC1  +NIIKPT  S  (  )  ) 
i  F  ( M  .01.  80)  f,0  TO  990 

(  Al  I.  DA  lAJ  I  TYP,N,M,/FP,RhMFP,TF>F  r  A,FR,VO|  |  N(.  ) 

IF!  I  TYP  .NF.  3  .AND.  I  S  I  GN  .FO.  MIN)  VOL  1  N(.  =-VOL  )  iv( 

II  YO=vo+vni i Mr 
WR I TF { 4,5? ) VO 

FORMAT!  '  OrilMPOTFO  R  F  SIJL  T  S  5  '  /  '  *  »  5X  »  '  VO|  i  IMF  '•  4  X, 
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00  2  N=  1  ,M 
! MOX ( 1 >=N 
AN=ARC<N1 
TN'  =  RHf)(N«R) 

on  ^  l=n.m 

I F { L  .FU.  00  TO  R2 
Al.  =  ARC  (  L  1 
TL=RHnn.»R) 

|  F(KFYMll-l  )  1  1  ‘tf  104  «  109 
lift  AHlllN.I.  )=0.n 
AM]  l(l.,N)=n.O 

]09  IF(KFYP11-1  M  n  .  11 ?tl  1 2 

HI  APll{MtU=0.0 

*  "  1  1  <  L  «  N  )  =  0  .  O 

t  ^(N,L)=0.0 
AP33(L»M>=0.0 

112  !N0X(2)=L 


92 


013630-9-T 


i if  i  103  .1  =  1  .  1  .  ? 

.H'A-  .1+  A 

•  -n  i  o*«  u  =i  t? 
i  =  3-u 

i i=inox<  u. ) 

I?=lNnx( i ) 

1  1  ?  =  RHfl (  I?,.IPM 

C\i  L  SP  Tl.PUFYPl  1  *  I  1  .  1?»JPA,  AP11  1  ,  API  3  3.  A- 1  ll.l) 

I  !■  i  K  F  Y  P 1  1  —  1  )  105, 10ft,  1  0* 

10-  AP11(  ll,]2)sAP11  (  11  ,  IPH-APil  1*11? 

A P 3 3  (  I  1  ,  l?)=AP33<  1  1  ,1?  1+APlV^T  \? 
i  f.  |F(KCVMll-r,107t10A,10A 

1  7  Ai.in(n,I?)sAMllCIl,I?)-AMIll4TI7 

!-ia  rriNTIMHF 
3  C  <  IN  T  I  Ml  IF 

CALL  SF  THP(kFypU  ,KFYM1  1  tNfi.,H,APl  1  1  ,AP]  AH,  AMI  \  \  ,  i  ) 
1MKFYP1  1  —  1  )10Ht? 09 ,704 
,  fir!  11  =  WO*  AP  II  1 

API  1  (N,L)  =  M  *  ( W  1  *  A  P 1 1  (  M  ,  I.  )+il*T|  ) 

API!  I  1.  »N )  =  AN*  ( W  1  *AP  1  1  (I.  ,N)+ll*TN) 

I I  =  w 0 <" A P I  33 

AP33(NfL)=AL*(Wl$AP33<M*!.  )+0*Tl  ) 

AP33< I  ,N)=AN*(W1*AP33<L.NJ-M)#TM 
?09  I FIKhYMl 1-1 ) ?] 0,3*3 

?\  O  A  M 1 1  I  l>i  t  L  )  =  A  L  •  <W1*AM1HM,L  J-AMI  11*W0=|  L  ) 

CALL  SFTIIPJ  1  ,0,L.N,H,API  11  tAP!33,AM(  11,0) 

AM  1 1  {  L,N)  =  AN*(W1*AMH  ( i.  ♦  N  )-AMJ  1 1  *WO*TM  ) 

OH  70  3 

«?  IF (KFYP11-1 )P3,fl4,ft4 

a  3  .i=ALOC-n  A.*TM/AN» 

A  P  1  1  (  A- »  M  )  =  ( 1 1—  1  »  0  )  *  AN 
A  P  3  3 1  f  *  »  N )  =  ( 1 1  •*- 1  #0)*AM 
PA  I F ( KFYM11-1 ) HSv3, 3 

p«s  ]F(F3  .FQ.  1.0)  C,n  TO  3 

no  fit-,  i  =  i ,  b 

»A  CALI.  SF  IMP  (  l  ,0,w,N,  1  ,U,n,ST3  (  !  )  ,  U 
0--.  5*1  1 .0-F=  )  *AN 

A»M  1  (  M,N)  =  P1  -II*  (WO*  {PHOUJ,?  )  *S1  3  (  ?)  +PH|,(  N,  5  )  3(  -  )> 

3  +W1  *(  PHO(  N,  l  )*ST3(  1  )  4PHfl(  N,3  )*ST3(  3)  +PHIMM,  A  I  I  j  l  4  )  + 

4  s>MO*  N,M  *ST3(  A  )  )  ) 

3  CON  TIN!  IF 

?  COM  INI  IF 
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II!  I  70  1  =  1  *M 

70  H(  !  )=TWnPl*RHII(  I  ,H) 

!F{KFYP1 1-1 171 

?.\  call  ofcompiaph 

CALL  SOLVFl API  1 ,X,H,M) 
no  77  1=1, M 

7?  X< I )=RHO(  l,P)*s?*X( I  ) 

CALL  I N TFG (  X  »NS6CS  , NIIMPT  S  » P 1  1  ) 

pi l  =  p 1 1 *pi /vo 

?4  1F(KFYMU-1 

75  CALL  HFCOMPl AK1 1 ,Mi 

CALL  SHLVFt  AP'l  1 
mi  7  4  1  =  1  ,M 

74  X(  n  =  RHIH  I  I  ,H)*XU) 

CAI  L  INTFG( X,NSFOS,NUMPI S,MJ  1  ) 
MU=P1*MU/Vft 

PH  I  f  (  KrYPl  L-l  137,45,45 

3?  Iin  79  1  =  1  ,  M 

79  H(  I  }  =  TW.  M*7I|  ,H) 

CAI L  OFCPMPl AP33,M1 
CALI  SOLVF( AP33,X,h,m) 
nn  134  I  =  1  ♦  M 

X  (  1  )  =Z  (  I  ,  R)*«M)(  I  ,  R  )  v  X  (  I  ) 

134  k  (  1  1  =  TWOP  1 

CALL  INTfcUl X,NSFCStNUMPTs,u33 ) 
CALL  S  01 V  F  (  A  P  3  3  ,  T  5  »  P ,  K- 1 
DU  35  1=1, M 
X (  I  l=RHO<  I tP)vT3(  1  ) 

35  H( I  ) =7 ( I ,8) -X (  1) 

CALL  lNTFr,(  X  .NSFPS  , MUMP TS, CAP  I 
r  a  p=  twup  i  ~c  a  p 

CALL  INTFM  P,MSFfA  ,A'i)MPi  s  »n  AM  ) 
r-AW=-TWnPI  4CAM/CAP 
P  3  3  =  (  T  W  n  P  1  *  P  3  3  -  C  A  P  "•  G  A  P  =  ( ;  A  M  )  /  V  (1 

lP(MRUn  -  1  )  39,39,44 

54  DP  35  1-1 tM 

1 F ( 1 -NF 1 ) 305,305,30 4 
305  M( |)= TWOP I 
(.1)  111  3  4 
304  m(  l  )  =  0,0 
3  5  ((lMTlNUt 
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(  A1  I  ‘sill  VF(AP33,T  1 
rill  307  1  =  1  fM 

30  /  X(l)  =  XHiiU,Al#Tlin 

f  A I  L  iNTbGJX.MsFG^fMUMPT*,,  f|.  ) 
r  A!  |  1  f"TF G(  X  ,NS1  .NIIMPT  S  ,TM  ) 

r.n  30R  1  =  1  ,M 

30«  X (  ! )  =  7  (1 ,H)*X(  1> 

r.A(  i.  !Mit-r;<x,N<^Gs,NUMPTs.tn 

!)H  TAP  =  -Ct  WOP  I  /VO)*  ( «l+OAM*Tl.  )**?/<  1  W-TwnP  i  *Tl  =  1 1  /C  aP  ) 

1 1=  P33+UF  l  T  AP 

30  wR  1  TF  (  ft,  AO  )  CAP,  GAM,  PI  1  ,P33 

AO  FORMAT!*  • ,SX, *CAPACI TANCE  s',HO,S/<  •,  5X  GAMMA  ‘  ,  7  X  ,'  =  ',  M  0.  S/ 

R  •  •  ,5X,'  Pll/V*  *7X,  •  =  •  ,F10.5/'  '  »5X, 'PH3/V' ,7X,  •  =  « 

I  F  (  Ni ROD  .FO.  ? ) D  R I T  b ( ft «  309  ) HE  L  T  A  P ,  O 

309  FORMAT! »  •  , 5X  » •  DEL  T  P  3  3  /  V  •  ,  ?  X ,  *  =  •  »  F  1  0  .  S  /  *  ',SX,  'rlS.IMI  P33/V=', 

ft  Fjn.ftl 

AR  JHKFYM1  1-1  )A?,  337,337 

A?  I  Th  (  ft,  A3  )  Ml  1 

A3  FORMAT!'  • ,5X,  "'ll /V' ,7X, »  = • tFlO,5) 

3  3  7  IFIIPPPT  « b  If  •  1  .AMI),  KFYP11  .EO.  0)  W«ITb(ft,<*M 
ft  (7(1 , «) ,RHO( 1 ,«) ,T3( 1) ,1=1 ,M) 

u  u  FORMAT!  1  0*  ,EX,  •  7  •  ,  10X  ,  'RHfl*  ,  1  2X,  •  T3'/ (  '  '  .  3  <  F  1  ?.  ft  ,  «;  X  )  )  ) 

r,n  to  37 

coo  '«'R  I  T*-  (  ft, ^9  1  ) 

09  1  FORMAT! »0*v*  ERROR  I  f';  DATA*) 

9Q9  CALL  SYSTEM 
FM(' 
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A, 3  SUBROUTINE  DATA  (IN,  MX,  MY,  ZEP,  RHOEP,  THETA,  B,  VOL) 

This  subroutine  is  called  once  for  each  segment  of  the  profile.  From 
the  input  specification  for  die  segment,  DATA  computes  the  ( z ,  p )  coordi¬ 
nates  of  the  necessary  sampling  points  on  the  profile,  the  quantities  cos  o 
and  sin  a  at  these  points  and  the  incremental  volume  of  the  segment. 

Arguments : 


IN 

Type  key  for  segment. 

MX 

Total  number  of  cells  in  segments  to  the  left. 

MY 

MX  +  ( number  of  cells  in  this  segment ). 

ZEP 

z -  coordinate  end  points  of  segment :  ZEP  ( 1 )»  z 
ZEP (2 ) «  z 2  . 

RHOEP 

p  -  coordinate  end  points  of  segment : 

RHOEP ( 1 )  *  p,  ,  RHOEP (2)  *  p  . 

1  2 

THETA 

Angle  ( in  radians )  subtended  by  a  circular 

arc  at  its  center. 

B 

Fractional  exclusion,  0  . 

VOL 

Incremental  volume  of  segment. 

Comments : 

Stored  in  COMMON  are  the  arrays  RHO  ( 80,  9),  Z(80,9),  ARC(80), 

C  ( 80 ,  9 )  and  S  ( 80 ,  9 )  which  contain  die  numbers  computed  by  DATA . 

For  the  Ith  cell,  the  subscripts  (I,  J)  correspond  to  the  points  s{ ^ 

of  ( 229 )  when  1  <  J  <  6 .  For  J  ■  7,  8,  9,  the  subscripts  ( I ,  J ) 

refer  to  the  points  s ,  ,  8,  and  s , ,  respectively  of  ( 211 ) . 

i  -  i  i+ 
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190- 

190? 

190*> 


SORROOTINE  DATA! I N, MX * MY , Z EP , RHOEP ♦ THET A, B, VQL ) 

DIMENSION  Z EP ( 2 1 , RHOEP ( ? 1 

COMMON  RH0(B0,9) ,Z(80,9) , ARC (BO) , C ( BO ,9 ) , S ( BO , 9 1 
DA  FA  STEP/. 3872988/ 

MXPi  =M,'  +  J 

FM  =  H.OA  T  (  MY-MX  ) 

IP(R  . NE .  1)  SOhSTP=.S-( 1.0-R  )«STEm 
1 F ( IN-? )1  ,2,3 
CO- 1.0 
00  TO  10 
CC=1.0 

S 1 ?  =  S  ?  N( fMFTA/2.0) 

A  =  /FP (? )-?  FP( 1 ) 

KAD=0.5-S0RT ( ( RHOF  P ( 1 ) -RHOEP ( ? ) ) **?  +  A*A  )/ST2 
On=A/ARS ( A ) 

T  =  CC»nf>*COS(  THETA/2.  )/ST? 

ZCNT=0.5v(ZbP ( 1 )fZEP( 2 )  +  T*(RHOFP( 1 ) -RHOEP (2  )  ) ) 

RHOCN  T=0. 5- ( RHOE  P ( 1 ) +RHOEP (2 ) +T* A ) 

02  =  7.  EP  (  2  l-ZCNT 
01 =Z  FP ( 1 ) -ZCNT 

VOL  =  3.  141593-ABS  <  A*  (  khOCNT**2+R AD*R  AD- «  02=?=*?  +01  *02+01  **?  )  /3.0) 

3  -CC*RH0CNT*(O?*(RH0FP( 2  J-RHOCNT )-0l* (RHOEP  ( 1  )-RHUCNT  )  +  k  AO'-R  AD 

3  *00$ THETA)) 

HE  TA=CC*f>OsTHFTA/FN 

1HF T1=ATAN2( RHOEP < 1 ) -RHOCNT , Z EP ( 1 )-?CNT ) 
n=ABc, ( KFTA^RAD ) 

K3=STFP*hETA 

Of i  90?  I  =MXP1  .MY 

PH  I  =  The  Tl  +  < l-MX-.SJ-K^TA 

IF(a  .FO.  1,0)  GO  TO  1 905 

00  1902  .1=1,2 

ANG  =  P“l+.5*U-1.5)*BETA*(  1  ,0  +  R  ) 

DO  1903  L=1 .3 
PSI=ANG+(L-2 l-SUBSTP^RFTA 
v<=  I.  +  3-  (  J-l  ) 

C ( 1 ,M)*-COSIN(PSl ) 

S { I ,M)=CC*CnS(PS! ) 

7 (I  ,M)=7CNT+RAn»CC»S( I .M) 

RHO«  I  , M )  =RHOCNT -CC*R  AD=>C  (  I  ,M) 

CONTI  MOF 
00  903  .1  =  7,9 
ANCsPHl  (  J-filvf-3 
C(  1  ,.))  =-CC*S!  N(  ANG) 

S  (  I  ».l  )  =CC*COS  (  ANG  ) 

7  (  I  ,.» )  =7f NT +R  AO*CC*S(  I  ,  J  > 

RHO(  I  .  J  )=R)inf.NT— CC*RAO*C(  !  ,J) 

ARC ( l ) =0 
R  F  TOR  N 


! 


013630-9-T 


3  DX=UtP<2)-ZEP(in/EN 

0Y=(RH0E?<2  )-RHOtPU  ) )  /EN 
t)  =S0R1 <DX*OX+DY*DY ) 

SI=DY/U 
C !  =  OX/U 

nn  917  I =MX PI* MY 
PH  t  =  FLOAT { 1 -MX )  —  .  5 
I F (  B  .EO.  1.0)  GO  TO  1800 
00  1802  J=1 » 2 

ANG=PHl4.5*( J-1.5)*<  1.04.,) 
nn  1803  i=i,3 
M=L+3*( J-l > 

PS! =ANG+(L-2 )*SUBSTP 

Z ( 1 ,M)=ZEP< l)+PSI*nX 
RHO ( I ,M)=RHOEP( 1 )4PSI*DY 
S ( I  ,M)=SI 
1803  C ( I »  M ) =C I 
180?  CONTINUE 
1800  00  913  J=7 ,9 

ANG=PHl4( J-B)*STEP 
Z(I,J)=ZEP(1)+ANG«DX 
RHO ( I  ,.J  )=RHOEP(  1  )+ANG*OY 
C(I,J)=CI 
913  S{I,J)=SI 
Q17  ARC(I)=U 

VOL=l. 047198* (ZEP< 2 )-ZEP<  1H*(RH0EP<1  >**2+RH0EP<  l!*RH0fcP(2)  + 
8  RHf)EP(2)**2) 

RETURN 

END 
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INTk  nunerlcally  integrates  quadratic  interpolating  polynomials 
approximating  the  data  on  each  segment  of  the  profL„.  When  the  profile 
is  composed  of  several  segments,  no  interpolation  is  performed  across  seg¬ 
ment  boundaries.  Hence,  the  integration  is  accural  even  for  disconnected 
segments,  e.  the  circular  arcs  of  two  spheres. 

\  rgu.su.:  its : 


V  Real  vector  of  function  values,  ordered  *s 

the  cells. 

NSEC  Total  number  of  segments  in  the  profile. 

NUMPTS  Integer  array  containing  in  NUMPTS  (I)  the 
number  of  cells  on  the  Ith  segment :  1*1, 
NSEG. 

SUM  Integral  of  V  across  the  profile. 


Comments : 

Stored  in  COMMON  are  the  arc  lengths  ARC  (I),  1*1,....,  N  required 
to  compute  the  integral. 
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«,l!«Rnt'TINE  I  NT  fir,  (  V.NSFG.NUMPTS  ,SUM) 

COMMON  PHDi 80,9 ) ,7 (80,9 ) »ARf (80 ) ,C<80,Q ) , 5(80, Q) 

niMFNSION  V  (  80 ,N(IMPTS  (  15  ) 

soM=o.n 

JAf.r.al 

no  8000  I  =  1 «NSPG 
T=A«C  (  .1  ACC  » 

I  =M(IMPT<;  {  I  ) 

Ms  I  +.1  ACC-1 

Sl)M,=  Sll,vl+T*(0.  6?5-(  V(  JACC  )+V(N  )  )-.  l?8-(  V(  JACC+!  )+V(W-l  )  I  ) 

I  F  (  L  /  ?  • N  F  *  (1+1)/?)  GO  TO  8001 

SliMrSUM+TiMO.hn  *>887*V(N-l  )-0. 0 8333333- V(N-?) +0.4  1oO‘,a/.  vt  )  ! 
8001  |M1=N-1 

Jl  Os  JACC  +  1 

00  800?  .l  =  .M.n,L«l  ,? 

800?  SIIK=SUM  +  0.3333333^T*(V(J-1  ) +4. 0*V  ( .1 ) +V  ( .1- 1  )) 

8000  ,|  ACC  =  .)ACC+L 

Pfi  TURN 
FNt) 
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Used  together,  DECOMP  and  SOLVE  solve  the  linear  system  AX  *  B 
DECOMP  performs  a  L  -  U  decomposition  of  the  N  >;  N  matrix  A  and 
SOLVE  performs  back -substitution.  These  routines  ars  adapted  from 
Forsythe  and  Moler  ( 1967,  pp.  68  -  69 ) . 


SUBROUTINE  OECOMp(UI  ,N) 

I  MENS  I  ON  UL(flO.fiO) 

COMMON  /SOL/IPSiaOi 
OH  6  1=1 *N 

ips<  n  =  i 

rv!Ml  =  N-l 

00  1ft  K  =  1 t NM 1 

RIE=0.0 

no  11  I  =K  »N 

I  P= I  PS (  I  1 

IF(ABS(ULUPtK)l  .16.  RIF)  00  TO  11 
RIF=ARS(0L( IP,K ) ) 

IOXPIV=I 

CONTINUE 

1MI0XPIV  .EQ.  K)  GO  TO  16 
.1=  I  PS  (  K  ) 

1PS(K1=IPS< IDXPIV) 

I P  S  <  I0XPIV)=J 
KP=IPS(K) 

P I  VO  T=U!_  <KP»K  ) 

KP l =  K  +  1 

no  ift  i=kpi,n 
IP=IPS( I ) 

BM  =  -lil  (  I  P,K  )/P  JVOT 
U|  (  I  P  » K  1  =  EM 
no  1ft  J=KP1 ,N 

Ul ( IP ♦ J )=UL ( IP.J)+EM*UL IKP.J) 
CONTINUE 
RETURN 
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SURROtlTINESOLVE  U)L  ,XtR,N» 
DIMENSION  UL(RO,RO),B<RO),X(RO) 
COMMON  /SOL/ I  PS (RO) 

NP l =N+ 1 

ip=ipsm 
x  c  1  >  =r ; ip) 
on  ?  1=2, N 
ip=ips<I) 

I Ml  =  I  — 1 
SUM=0. 

no  l  j= l ♦ i mi 

SUMs  SIJM+UL  UP(J)*XI  Jl 
X ( I >  =  B( I P  >  — SUM 
I  P=  I  PS  (  N I 

X(MI  =  XIN)/ULUP,N) 

DO  4  IBACK=2tN 
I=NPl-IRACK 
IPs  I  PS (I ) 

I  PI  =  I +1 
SUM=0.0 
on  3  JsIPUN 
SUM=SUM+UL( IP«J)*X(J) 

X(  I  )s(X(  n-SUM)/(JL  (  IP*I  ) 

RETURN 

END 
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A. 6  SUBROUTINE  ELLI  (Ml.  K,  E,  KPR.  KEY) 

This  computes  the  elliptic  intsgrals  K(m)  and  E(m )  and  the 
derivative  K'(m)  from  their  power  series  approximations  ( see  Section 
4.2). 

Arguments : 

M 1  Real,  the  quantity  (1  -  m ) . 

K  Real,  K(m). 

E  Real,  E  ( m ) . 

KPR  Real,  X'(m). 

KEY  Integer:  0  Compute  K,  E  and  KPR; 

1  Compute  K,  E  but  omit  KPR . 


SURRDUTINF  ELI  H Ml » K , E ,KPR , KEY  1 
PEAL  Ml,KtKPR 
T=-Ainr,iMi ) 

K  =  \  .3R*?94+.5*T+M1  #(9  .*F>f>344E-?+.l?49R59*T  +  Ml  =  (  3. 59009? E~? 

5  ♦  6.  ftfiO?4q  £-;>*T+Ml s  (  3»  7425ft4E-?  +  3«.3?R35<>E-?:i‘T  +  Ml »  ( l  •  451 1  9fcF  -? 

5  +4.4l7R7fi-3#T  n  n 

Ml  ,0+Ml*  l  .443751 4+. 2499 837*? +Ml*<  ft.  2f«0F>01  E-2+9. 200 1 *F-2* i  +  M1*{ 

R  4.7573H4E-?+4.0ft9ft9  8F-?*T«-Ml*  (  1  .?3ft60ftE-?+5.?ft449ftE-3*T  11)1 
IF {KEY  .FO.  1)  RETURN 

KPR=.S/M1  *  ?.R3??‘jf-?  -  .1249R59M  *  Ml  * ( -? ,  999362F-3-.  1  3  7*05* T 
ft  +Mloi-7.fi9933ftF-?  -  9  .9R50ftftF-?*T  +  Ml  *  ( -5.  3ft299Bfc-?  - 
5  1  .76714RF-?  *  T  1)1 
RF  TURN 
FNU 
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A. 7  SUBROUTINE  SETUP  (KEYP11,  KEYM11,  I,  J,  L,  API11, 

API33,  AMI  11,  IJ) 

This  is  essential  in  computing  the  linear  systems.  Specifically 
SETUP,  after  calling  ELLI,  computes  the  quantities  API  11  (K1  of 
eq.  182),  API33  (K0  ofeq.  178),  (eq.  196)  and  (eq.  197). 
The  quantities  Q  and  Cl  are  used  to  compute  f(i,  j)  (AMI  11)  of 

A  C 2 

eq.  (227). 

Arguments : 

KEYP11  0  when  computing  API  11  and  API  33, 

else  1. 

KEYM11  0  when  computing  AMI  11 ,  else  1 . 

I  Subscript  of  observer  (unprimed)  cell. 

J  Subscript  of  remote  ( pruned )  cell. 

L  Index  of  the  point  within  remote  cell  for 

which  the  kernels  are  to  be  computed 
( see  DATA ,  Comments) . 

APIil,  API33  Described  above 
AMI  11 

IJ  0:  use  last  value  of  Ml  in  kernel 

computations ; 

1 :  compute  new  Ml . 
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SUBROUTINE  SETUP (KEYP!1,KFYM1 1,1  ,  J  , L  ,  AP  i  1 1  ,  Af>  1  33 ,  AMl  H  .  1 .1) 
COMMON  ft  HO  (80 .9  1,7.  (80,9)  ,ARC  ( 80 ) , C ( RO ,9 ) , $ ( 80, 9 ) 

PEAL  MfMl,K,KPR 
70=Z(J,L )-Z( !,R> 

R=RHO( 1,8} 

RP-RHO< J,L ) 

I  F ( I J  .FO.  0)  GO  TO  115 

RRP=R*RP 

A 1 =RR  P+RR  P 

A2  =  R<‘R  +  RP*RP  +  Zn*Z0 

Ml=(A?-Al)/( A2+A1  ) 

M=1 .-Ml 

CALL  FLL1 (Ml ,K,E,KPft,K£YMlli 

AO=M/RRP 

A1=S0RT( AO) 

A?=M+M 
A3=  2  • -M 

!F(KEYPli-l  )in, 114,114 

113  API ll=Al*(A3*K-E-E) /M 
API33=A1*K 

114  A  I s  A 1  *  A  0 
A3=.5*A3 

1F1KFYM11-1 1115,116,114 

115  A  0  =  C  ( .) »  L  ) 

0^1  =— A1  * (  ,25*K-A3*K.PR  ) 

0^2=A1* ( F-A3* ( ( 63+M)*K  -  A2*A3*KPR ) )/ (M*M) 

AMni  =  ft*AO*OM?  +  (ZD*$(  J,L)-RP*AO)*OM1 
114  RETURN 
FNO 
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